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CHAPTER I 


Introduction 


One of the primary goals in this dissertation is concerned with the develoi)nu'iit 
of robust hybrid finite element-boundary integral (FE-Bl) techniques for modeling 
and design of conformal antennas of arbitrary shape. Both the finite element and 
integral equation methods will be first overviewed in this chapter with an emi)hasis 
on recentlv developed hybrid FE-BI methodologies for antennas, microwave and 
millimeter wave applications. Xhe structure of the dissertation is then outlined. We 
conclude the chapter with discussions of certain fundamental concepts and methods 
in electromagnetics, which are important to this study. 

1.1 Overview 

The development of simulation techniques for conformal antennas typically mounted 
on vehicles is a challenging task. By and large, existing analysis and design meth- 
ods are restricted to planar and mostly rectangular patch antennas. These tech- 
niques have difficulty in being extended to non-rectangular/non-planar configurations 
loaded with dielectrics and comprised of intricate shapes to attain larger bandwidth 
and gain performance [1—4]. Moreover, practical antenna designs ma\ also ie(|Uiie 
a sophisticated feeding structure, such as coaxial cable, microstrip line, stripline. 


1 



|>ro\iinity or a|>crturc rou])h‘(l cirrnit ruMwork. rw . <uid rcjiuition niriiuKU 

are not easily adapt ahl(‘ to inodclino i1k‘s(' st nu t uivs. (‘>p(‘t i«dlv in iIk* |)ivvriu (' 
of finitely sized diedectric loadine-. Partial diffcMeiit ial <*(juation iPDI') t <'cliiii(|U<> 
(e.g. finite element and finite difference methods) ma\ al>o (\\peri(Mic(‘ difiic ultieN in 
modeling unbounded field problems, such as thos(' found in antenna radiation and 
scattering. The motivation of this dissertation is therefore based on the iu'(h 1 to (h*- 
velop general-purpose analysis techniques which can accurately simulate conformal 
antennas of arbitrary sliape with diverse feeding schemes. With the rapid growth of 
personal cellular GPS and other communication systems, there is an increasing need 
for such techniques since even traditional and protruding low frequency antennas 
are being re-designed for conformality and to meet requirements for a host of new 
applications [5,6]. Besides, most development of computational electromagnetics in 
this subject can be applied to medical diagnosis and treatment which have shown a 
tremendous research and application potential [7,8]. 

The complexity of new antennas demands that analysis and design software be de- 
veloped based on methodologies that are robust, versatile, and geometrically adapt- 
able. Recently it has been demonstrated that the finite element method when cou- 
pled with the more traditional integral equation approach becomes quite attractive 
for modeling a wide variety of existing and emerging antenna configurations [9]. The 
finite element method is indeed ideal for modeling the interior volume of the an- 
tenna structure (multi-layer substrate, finite size dielectric loading, stacked element 
design, feed network and cavity volume, etc.) and is one of the most celebrated 
analysis methods in engineering. On the other hand, the boundary integral offers 
the most accurate representation of the fields exterior to the antenna. Thus the com- 
bination of the finite element and the boundary integral (FE-BI) methods provides 



for tlic liandline, (jf the o(>onu‘t rical coin|)l»‘xit> witlKnit ci.)in|)roinisni!: ;i< cmary. 1 hi- 
livbrid ni(‘lhodology appears to be very attrartive for coiifonual antenna mod( liny. 
However, its development and apjdication to inor(‘ ])ractiral and enu'rginy antenna' 
presents us with many theoretical and numerical challenges, which will \>v ('Nlensi\ cly 
investigated in the work. 

Specifically, mesh termination plays an important role in F L.M simulations and. 
in manv cases, the accuracy is subject to the performance of tlie domain truncation 
scheme. For conformal antenna modeling, a boundary integral (Bl) equation has 
been employed in this dissertation for terminating the antenna's radiating surface 
and this method is theoretically free of approximation. Thus, a desired accuracy can 
be achieved without fundamental limitations. Antenna configurations of arbitrary 
shape can be readily tessellated using mesh generation packages in the context of 
the FE-BI technique. In modeling the interior region or the feed network, a superior 
artificial absorbing material — perfectly matched layer (PML) — has been used to 
ensure a minimum impact due to truncation walls. An intensive study of the PML s 
performance has been carried out and the optimal selection of PML parameters has 
been designed and employed herewith in shielded structure modeling. 

Frequency domain methods provide the necessary information for engineering 
design. However, when wideband responses are needed, they can quickly become 
expensive compared to time domain techniques. A method, referred to as the asymp- 
totic waveform evaluation (AWE), can be used to alleviate this issue. It has already 
been successfully used in VLSI and circuit analysis. In the context of the FEM. we 
shall investigate the suitability and validity of AWE for simulating MMIC devices. 

One of the important issues in antenna analysis is the feed design. Modeling a 
feed using the finite element method is indeed a challenging problem, and a sim- 



[>lili('ci i>rol»e feed iiiodcl fail- lo accurately [)redict the iiipui inii>edaii( i'. On iln 
otlier iiand. tlie numerical system can l)ecome ill-cumlil ioneil uIkmi a feed neiwuik 
is modeled without careful considiTations. In tlii> dissertat it)ii various feed iiu.drU 
will be investigated in consideration of accuracy and efficiency. They inclmh' c uiK'nl 
and voltage gap generators, stripline, microstrip line, coaxial cable, aperture cuu|)led 
microstrip, etc. 

In regards to the development and applications of the simulation technic|ues. the' 
test and design benchmark models of particular interest are microstrip (rectangular 
and circular) patch antennas, dual-stacked patch antenna, ring slot antenna, and 
cone antenna, etc. It is noted that some of them are not necessarily planar or 
conformal. 

Referring to the dissertation structure, we begin with a description of electro- 
magnetic fundamentals and then proceed to discuss the boundary conditions, equiv- 
alence principle, Dyadic Green’s functions and the related theorems. The finite 
element method as applied to time-harmonic electromagnetic fields and waves is 
subsequently described and the basic FEM equations are derived from both varia- 
tional and Galerkin techniques. The derivation is given in algebraic form allowing 
the inclusion of general anisotropy. The emphasis of the discussion is on the gen- 
eralization of the variational functional and Galerkin techniques when anisotropic 
and lossy materials are present. Chapter 3,4,5 and 6 discuss the development of 
edge-based FE-BI techniques with significant efficiency improvement for antennas 
and feed network modeling. The emphasis in these chapters is on developing novel 
methodologies to minimize the required computing resources. 

Chapter 7 is devoted to circuit modeling where specialized truncations suited 
for guide wave structures are presented. The perfectly matched layer (PML). an 



aiiisot r(i[>ic artifinal absorber usc'd for iiu'sli Iruiicat urn. i> nufsl iiiated in leini' ul 
j^erforinaiice ami a|)i)licalioiis. 

\\ icleband svstem responses j)iompt us to look at more c'llu lent analysis tooU 
to replace tlie current brute force frequency domain analysis approaches. ( haptei S 
discusses a preliminary development of the FEM in connection with the .\\\ L. 

In the last chapter, we summarize and discuss the anticipated future research 
work to extend the capability and applications of the robust FEM develoi>menI , .A 
list of suggested topics is included with specific recommendations. 

1.2 Fundamentals of Electromagnetic Theory 

Since many fundamental concepts and theorems of electromagnetics will be em- 
ployed, we will describe the pertinent ones in this section for reference purpose. 
This will also ensure consistency in nomenclature and conventions throughout the 
dissertation. 

The vector wave equation — the only partial differential equation (PDE) con- 
sidered in this research — will be first derived from Maxwell equations. Various 
boundary conditions will be studied to establish the general mathematical models of 
boundary value problems (BVP). The equivalence principle, uniqueness theorem and 
the half-space dyadic Green’s functions are then briefly discussed for EM solutions 
in radiation and scattering problems. 



1.2.1 Maxwell Equations 


1 iiiK* harmonif Maxwell ccjiial ions of diffen'ntial form in a linear, ani^oi ioj>n ,iml 
uniform medium are given by [10] 


V X E 

V X H 
Vf -E 
V H 


• H — M, 
j^t ' E “f" «J I 

Pe 

pm 


.:i) 

.4) 


where E and H are the electric and magnetic field intensity, respectively. is the 
radian frequency and the factor is assumed and suppressed throughout this 
dissertation; M, and J, are the impressed magnetic and electric current, respectively, 
to serve as possible sources in the medium under consideration; finally and p„, 
denote the electric and magnetic charge density. Both M, and pm are fictitious and 
non-physical quantities, which facilitate the formulation of physical problems when 
the equivalence principle is employed. The material tensors t and p represent the 
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with to and po being the free space permittivity and permeability. 



Tlie procedure to derive the vector wave ec|uation l>e|iiiis 1>\ eliniiiial iiiu our ot 
the two field quantities from (1.1 I and (1.2). To do su, we lust take a dot prodmt 
of (1.1) with the tensor p”' and then lake the curl on hotli sitles i>f (1.1) to ohtam 

V X /7~' • V X E = X H — V X Pr ■ ( 1 ■ I ) 

Substitution of (1.2) into (1.7) yields 

V X ' • V X E^ = jJotolr ■ E — — V X 

or 

V X • V X e) - k^fr • E = - juJUoJ, - V X ■ M,) (1.8) 

where ko = is the free space w'ave number. The dual of (1.8) is given by 

V X ■ V X h) - • H = - jcJCoM, + V X (i;^ • J.) ( 1 .9) 

and can be similarly derived starting with (1.2). Equations (1.8) and (1.9) are the 
vector wave equations of the desired form. 

1.2.2 Boundary Conditions and Boundary Value Problems 

Three types of boundary conditions are typically encountered, and in the context 
of the finite element method, these boundary conditions must be considered and 
carefully treated. In what follows we shall discuss these conditions. 

Dirichlet Boundary Condition 

Consider two media separated by a surface T whose unit normal n points from 
1 to medium 2. The fields on two sides of the interface satisfv the relation 


f) X (E2 — El) = —Ms 


( 1 . 10 ) 



s 


whrrr M5 is rj firtitious magiielic surface* curre'iU and Ei and Ej arr ilir (‘ItH tiu tifld 
inside medium 1 and medium 2. respect ivedy. If nu’dium 1 is a p(‘if(‘clly maiiiiriK 
conductor (PMC), then Ei vanisli(‘> and (1. 10) 1k‘cuiiu*s ft \ Ej = — M^. I lu* ‘'urf.ur 
magnetic current can either be an impressed source ((‘xcitation ) or may r(*|)r<‘>eiit 
a secondary (induced) current. If medium 1 is a perfectl\ edectric conductor (Pl\('). 
fi X E2 also vanishes and thus M5 = 0 on the PEC surface. 

Similarly, for the magnetic field, 

n X (H 2 - Hi) - J, (1.11) 

where denotes an electric surface current. The PEC surface can support electric 
currents, given by n x H2 = J5, since Hi is zero within the conductor. By duality, 
the PMC surface does not support electric currents, i.e. V x H2. 

The relations (1.10) and (1.11) are inhomogeneous Dirichlet boundary conditions. 
They become homogeneous when M* = 0 and = 0, and in those cases they imply 
the tangential field continuity across the dielectric interfaces. Often, and are 
introduced as fictitious currents when applying the equivalence principle (except in 
special cases where they are specified a priori). The implication of this issue will be 
discussed later in the development. 

Neumann and Mixed Boundary Condition 

In formulating a physical problem using hybrid finite element methods, we usually 
work with either E or H field. If, for instance, w'e choose E as the working quantity, 
then (1.11) must be rewritten as 

n X V V X e) — (p V V X e) j = — ju;Js (1-12) 

where (p ^ • V x e) ? = 1,2 are evaluated just inside the ?'th medium approaching 
the boundary (from the ?th medium). If medium 1 is a PEC, then V x Ej = 0. and 




(l.rj) roclucos to a slandarci Neumann houiuiary rumlii i(.>ii. l*y whicii a luii-naiiit 
on tlir derivative of E at the interface is defined. I lie dual of i" ei\t'ii hy 

i, X f(r’ • V X h) ^- (r‘ -V X h)_] = 

and this condition is used when working witli the H field. In many applications, the 
single field formulation is often desired since the system size may be kc’pt minimum 
in this manner. However, it is already seen that the single field formulation iiiijilii's 
use of the second order conditions referred to as natural condition.-<. l ortunately. it 
is rather straightforward to impose this type of conditions in regard to finite element 
simulations. 

As for mixed boundary conditions, an example is the resistive surface where the 
electric and magnetic fields satisfy the condition 


fi X n X E + Rn x [H] 


0 


( 1 . 14 ) 


= H+ - H" the field 


with R being the effective resistivity of the surface and [H] 
difference above and below the surface. This is a typical mixed (third type) homo- 
geneous boundary condition. Another example of a mixed condition occurs in tians- 
mission line problem (e.g. a coax cable, or other guided wave structures), where the 
electric and magnetic fields at a cross-section of the line are given by 




E = E’e'^-’ + rE‘e 


H = 


( 1 . 15 ) 

( 1 . 16 ) 


and 


fi X E* = -ZH‘ 


( 1.17 


where n = -z and (E', H') are the incoming fields before encountering a discontinuity 


or load along the transmission line. .Also. Z is the wave impedance associated with 
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tlic* transmission lino mode of the frui<l<‘ \va\o >tni('tur(\ Kliminal inc 1 fiuin i 1 . 1 t 
and (I.IG) in \ ic\v of (1.17) yields the rolation 

// X E ZH = 2// V E' f ~ ‘ ' ( 1 . 1 > I 

which is an example of inhomogeneous mixed (third type) boundary condition. 
becomes apparent when H is expressed in terms of the derivative (curl) of E.^ and 
therefore the left hand side contains both differentiated and undifferenliated cjuan- 
tities. (In this case, the right hand side is usuallx' considered as a known function.) 
The mixed boundary condition (1.18) is found very useful when applying the FEM to 
guided wave structures for truncation and excitation simultaneously. It is basically 
a form of absorbing boundary condition (ABC). 

1.2.3 Uniqueness Theorem and Equivalence Principle 

The uniqueness theorem and the equivalence principle will be explicitly or implic- 
itly applied to this work when dealing with integral equations to terminate the FEM 
mesh and when evaluating the far-held pattern. Together with dyadic Green's func- 
tions, it becomes convenient to apply these concepts to construct integral equations 
associated with various geometries in radiation and scattering problems. It is our 
intent to discuss the theorem and the principle (without proof) for later applications. 
Uniqueness Theorem 

Partial differential equations (PDE) can be solved using various approaches and 
the corresponding results can also be represented in numerous forms given certain 
boundary conditions. Moreover, many (boundary, initial, natural, essential, radi- 
ation, etc.) conditions of PDE models can be extracted from the mathematical 

^Care must be taken when a curl operation is performed at a boundary discontinuity. It should 
be appropriate to evaluate the field derivative at a distance from a discontinuity and then let the 
distance tend to zero. 


s[>crifirations of well (lelincd ])liysical l)t■olll^•ln^. 1 ho <nu‘>lioii then ari"i-' ;i' u> liou 
to relate the solutions and how nian\ eonditions ar»' suflieieni tu achiew the lot- 


reel ' solution. I ni(|ueness t heorcins offer t he answer to this ([iiest u)ii. Sj>e( itic all \ in 
electromagnetics, tin EM f-ohttioiis an uniqinly ddtrmnnd by tin .•.oiiirt - ni a yivtn 
region pim the tangential components of the electric field on boundarns. or plus tin 
tangential components of tin magn et ic field on boundarn s. 


Equivalence Principle 

From the uniqueness theorem, an EM problem can be uniciuely solved if the 
tangential (either the electric or magnetic) field component at the boundary is jire- 
scribed. In this work, of interest is an EM problem where a dielectric inhomogeneous 


region exists in the presence of a large PEC platform, probably coated with a dielec- 
tric slab. The typical geometries are shown in fig. 1.1. where we consider the upper 
half space to be the exterior region and the cavity the interior region. 

In EM analysis, the fields in the exterior region can be represented in integral 
form containing the equivalent current sources. From (1.10) or (1.11) the tangen- 
tial electric or magnetic field near the aperture (or the discontinuity region) may be 
equivalently expressed in terms of the surface currents M, and/or J,. By ‘equiva- 
lence’, we demand the field distribution remain the same when the fictitious surface 
currents are used to replace the interior region (cavity volume). It can be shown 
through the uniqueness theorem that this substitution indeed ensures an identical 


EM field distribution in the exterior region. 

When the interior region is excluded from consideration, the current sources in 
(1.10) and (1.11) may be arbitrarily chosen leading to an infinite number of choices 
for the equivalent currents. However, in our work the field behavior in the intcrioi 


-See the proof in reference [10]. 



equivalent currents 


Figure 1. 



ground plane 

(a) 



cavity 


coated ground plane 


(b) 

: (a) Recessed cavity in a PEC ground plane, (b) Recessed cavity in a 
dielectrically coated PEC ground plane. 



l.i 


rccioii is also iK'cdcd. and most sju'dlically. tlir cuuidiiiii of llic ti<dds in I hr innn <ind 
outer cavity rrgions IS d(“si r(’d . It is t iKTcforc c (.>nv('nicnt tost'l<'cl tlic total l<inc.<'iit ml 
electric or magnetic field to specif} the ec|ui\aleiit currents as 

Mj = E X h and Jj = >i x H. ( l.lfM 

This choice implies the assumption of zero interior fields when the ('Xterior region 
is considered, and zero exterior fields when the interior region is needed, h ig. (l.d) 
and (1.3) illustrate the details of applying the principle, where the fictitious curix'iits 
affect the region of interest (ROI) only with zero EM fields outside of the ROl. It is 
observed that this choice of equivalent currents permits a convenient interior/exterior 
system coupling for the ‘Total field formulation in hybrid FEM applications. 

1.2.4 Integral Equation and Dyadic Green’s Function 

The Dyadic Green’s functions are particularly convenient for constructing integral 
equations in the presence of certain canonical platforms. For a planar structure, the 
platform of particular interest is the PEC infinite ground plane in which a cavity is 
recessed with dielectric loading or absorption depending on applications. 

The choice of the dyadic Green’s function varies depending on the FEM formu- 
lations. For the electric field formulation, we are seeking an appropriate integral 
representation to find the magnetic field in the exterior region using the information 
on or near the region of the aperture. To this end, let us start with the structure 
containing a possible protrusion as shown in fig. 1.4, where the equivalence princi- 
ple has been used on the outer contours of the structures to obtain the equivalent 
currents. 


Consider the wave equation 

V x V x G(r. r) - i^’VoeoG(r, r') = -I(5(r - r') 


( 1 . 20 ) 





ground plane 


Region of Interest (ROI): Exterior 

(a) 

equivalent currents 



Region of Interest (ROI): Interior 

(b) 

Figure 1.2: Illustrations of equivalence principle when applied to the structure shown 
in fig. 1.1a. 
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equivalent currents 




{E,H)=(0,0) 




ground plane 


Region of Interest (ROI): Exterior 

(a) 


equivalent currents 


cavity 



(E,H)=(0,0) 




Region of Interest (ROI): Interior 

<b) 


Figure 1.3: Illustrations of equivalence principle when applied to the structure shown 
in fig. 1.1b. 




(a) 



(b) 


1 . 21 ) 


Figure 1.4: Examples of protruding configurations (a) on a planar platform (b) on a 
curved platform in consideration of the equivalence principle. 

where G is the dyadic Green’s function G in association with (1.9) (assuming M. = 
0), and I is the idem factor defined as I = ix yy + zz. Also, note the identity 

III ■ (V X V X Q) - (V X V X P) • Q} dV 

= - IJ n • [P X V X Q + (y X P) X ^ dS 
and upon setting P = H and Q = G, we get 

III • (V X V X G) - (V X V X H) • G} dV 

X y X G + (V X H) X G] dS 

From (1.9) and (1.20), the left hand side (LHS) of (1.22) reduces to 
LHS = -H(r') G(r|r')dV 

and the right hand side can be rearranged as 


:i.22) 


RHS = / / H • [h X y X ^ + (y X H) • [h X G] dS 



(\:2-U 


K(|uatiiig tlic LHS and HUS yields 


H(r) 



J ■ G(r'jr) ' 

n X V' X Gj + ( V' V H ) • [u x Gj } <l>' 


where r and r' have been interchanged without loss of generality. .-\s can he realized. 
V is the volume containing llie distributed electric current source and 5' is tin- 
surface enclosing the entire upper half space. 

To eliminate the curl on J. we use the dyadic identity. 


V • (J x G) = (V X J) • G - J • (V X G) 


and the divergence theorem to get 

J . ( V' X G) dv - yy .■> . (j x g) </x ( i ..m > 

where the Sommerfeld radiation condition was invoked to eliminate the integral at 
infinity. Therefore, S' is only over the outer surface of the body. 

It remains to represent the surface integrals in terms of the electric field near the 
cavity since this field is typically the computable quantity. This is carried out b\ 
inserting (1.24) into (1.23), yielding 


H(r) 


~/// ^ G)dV 

- {h X V' X G) + (V' X H - J) • (n x G)} dS' 

~ III 

- JJ {H • (n X V' X G) + ju;eoE - (n x G)} dS' 


T.25) 


where the Maxwell equation (1.2) has been used. It should be remarked that the 
above field representation is general, i.e. not restrictive to planar or conformal cases. 



1^ 


I'or instaiirc. in ihv ])r<‘>cnr(‘ of a PKC‘ platfunn. i'- xalid fur ludiiiu (\m 

figurations as shown in fig. l.-l. lii tli<ss(‘ cas(‘>. tlu- suifarc iniuoral ioii^ an* (.init'd 
out o\er the platform plus the ouU‘r conlours of the strueturi's. 

The field representation (1.25) shall l)e examined and compar(‘d for (‘oidormal 
and protruding structures. To this end. we rewrite the surface integral as 


^ sur fact — 





(V X G) • (5 X H)dS' + X (V' X GV • (5 x E) } 




V 


IL {» 


X G)clS' - -^{i) X E) ■ (V' X G) }► (IS' 
u-y/o 


./lit 


where T denotes a transpose operation of the dyadic and the integral in the last 
step is proportional to the electric field. If G(r|r') is the electric dyadic Green's 
function of the first kind defined as n x G = 0, the first term in the integrand 
of (1.26) vanishes on the platform, provided S' is coincident with the platform. For 
dielectric protrusion, this term reduces to the integration only over the outer contour 
not conformal to the platform. An alternative is to define an electric dyadic Green's 
function which satisfies the condition h x (V' x G)^ = 0. As can be seen, this 
definition of the Green's function equivalently leads to the same vanishing term in 
(1.26). G is referred to as the dyadic Green’s function of first kind. The equivalence 
of both definitions can be proved from the symmetry properties of the dyadic Green's 
functions [11]. 

For a planar PEC platform, G reduces to 


G(r|r') = Go(r|r') - Go(r|r') + 2f ~G'o(r|r') 


(1.27) 


w'here Go is the free space Green’s function given by 

,-jA-o|r-r'i 


G'o(r|r') = 


47r|r — r' 



1 " 


and 


Go(rjr') = 


Inserting (1.2G) and (1.27) into (1.25). we obtain 

H(r) = H'"^'(r) + H^'^(r) + 2jkYo jj^ Go(rir') • (n x EldS' (1.2S) 


This is the desired form of the magnetic field representation used to estal)lish the 
boundary integral equation for a planar platform. 



CHAPTER II 


Finite Element Analysis in Electromagnetics 


The finite element method (FEM) has been applied to electromagnetics (EM) 
since several decades ago [12]. Especially in the late eighties and early nineties, 
it is observed that the publication volume associated with the FEM in electrical 
engineering grew in a fairly rapid pace [13]. This is primarily because electromag- 
netic problems in engineering designs become increasingly complex and analytical 
approaches or other numerical techniques no longer meet practical needs. With its 
numerous attractive features over other numerical techniques, the FEM has been 
extensively investigated and exploited for various EM applications [13]. 

This chapter is organized as follows. Section 1 and 2 describe the theoretical 
formulations to construct the FEM equations. These are usually considered the 
indispensable fundamentals of the technique, even though some interesting issues 
associated with these basics are still in development stage, especially in terms of nu- 
merical implementation. Of interest in this context is the discussion of the variational 
functional and Galerkin's techniques when applied to general anisotropic and lossy 
electromagnetic problems. This topic is one of the least studied and documented 
in the literature related to computational electromagnetics. .Anisotropic materials 
have been used for domain truncations (refer to Chapter 7) and therefore the general 
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v(“Clor or U'lisor form will !>«' used wlieiievfr i>i>s>il»l<‘ iii the iieees>ai\ deii\<ilioii' 
for tlie l EM. \\ ith this type of fornuilatioiis. isotropic system^ ina> l>e regarded a- 
special cas<*s. 

The chapter is concluded with the discussion of llie physical quantities tor antenna 
analysis in association with the computation of electromagnetic fields. 1 he formulas 
ffiven in this context require minimum amount of effort for computations. 


2.1 Functional Formulation 


The FEM was first developed with the aid of functional analysis. Tradition- 
ally, many standard boundary value problems (B\T) encountered in practice can be 
equivalently related to the extremization of a certain variational functional. W ith 
the Rayleigh-Ritz procedure to project a continuous function space onto a discrete 
finite expansion space, the variational functional method can be used to solve those 
physical problems and therefore becomes one of the two important approaches to 
formulate the FEM. A functional version of the FEM for the vector wave equations 
(1.8) or (1.9) is discussed in this section, which can readily incorporate boundary 
conditions, sources, resistive cards and other constraints into the formulation. It is 
regarded as a natural, convenient and sometimes physically meaningful approach. 
Furthermore, the functional may represent a true physical quantity (e.g. in low 
frequency, power transmission applications) and hence this formulation provides a 
feature of merit for its evaluation. Also, as can be seen in this chapter, the varia- 
tional method in general non-self-adjoint cases may be rigorously treated to result 
in a final symmetric system, a subset of which is identical to that obtained from 
Galerkin’s technique. Last, but not the least, the variational functional formulation 
can be used to validate the expressions based on Galerkin s method. 





( oiisidcr cl typicrti radiation oi scattt'niiu, i)rol)!('tii >lio\vn m liu. _’.l. uhrif ilir 



Figure 2.1: liiustration of a typical conformal antenna configuration. 

radiating elements (or array) are enclosed in a region $7. The platform surrounding 
the radiation/scattering geometry can be a planar ground plane, certain canonical 
shape (cylinder/sphere), or even a doubly curved surface in which case the Green's 
function is not available. In fl, electromagnetic fields satisfy the wave equation (1.8) 
or (1.9), which can be concisely described using a linear operator £ given by 

£^ = K, (2.1) 

where $ denotes the field E or H, and 


£ = 

(vx^;^vxj 

1 - (^o^r-) 

for electric field 

(2.2) 

£ - 

(v X !;' • vx) 

- (^O^r-) 

for magnetic field 

(2.3) 


K, is the source term associated with the impressed electric and magnetic currents 
and may be explicitly given by 

K, = — —V X for electric field ('7.4) 

K, = -f V X for magnetic field (2.o) 





As alrcadv incntioncd. £ is a linear operalor and <.>n(' ciui readil> slu>\\ iliat loi 
synimotric dielerlric tensors /7 and T. tlie pertinent iunetional of the original I’DI. 
lias the form 


JT((})) = - < 4>.£4> > - < ^.K, > 


CJ.(d 


where the inner product <. > is defined as 

<A.B>=fAB-d\ (£7) 

Jq 

(with B’ being the complex conjugate of B) for lossless media, or more generally as 

<A.B>= [ A BdV (2.S) 

Jq 

for both lossless and lossy media. 

The equivalent variational problem can now be stated as the extremization of the 
functional (2.6) in conjunction with the essential boundary conditions (e.g. Dirich- 
let EC’s). Specifically, the boundary value problem is equivalent to the following 
variational model 


= 0 


(2.9) 


Essential Boundary Conditions 

Because the effect of complex materials on the resultant system is of primary 
interest to us, we restrict most of our discussions in this chapter to homogeneous 
Dirichlet boundary’ conditions unless otherwise specified. Therefore, the variational 
approach of (2.9) ensures a symmetric numerical system. This is significant since 
many physical problems retain a certain symmetry property and the corresponding 
mathematical models should therefore reflect this property. Moreover, the symmetry 
of a numerical system is always desirable since it leads to more efficient solution and 
less storage requirements. 



I Ik- <'.\isu-iK'c of the fiinrtioiial CJ.(i) n‘<iuir«-> the operator € l>e x // mljniiit. ,t 
proiKTty usually defined to satisfy 

< £<I>. 'I' > = < <I>. £'!' > I d.ioi 

where <I> and 'J' represent any two admissible functions. If (LMO) holds, not onl\ 
does the numerical system derived from the functional (2.()) remain syminetrit . the 
minimization/maximization becomes physically meaningful. 

Of most concern is the situation where the partial differential operator of a system 
is no longer self-adjoint. Mathematically there exists no such a functional in 

the case similar to (2.6). A typical example is the presence of a lossy and anisotropic 
medium, whose dielectric material tensors are not symmetric or Hermitian. and this 
type of problems is more often seen nowadays. The development of finite element 
methods for those problems is still at an early stage because it involves numerous 
challenges. 

Traditionally, these physical problems were fictitiously simplified and dealt with 
using available numerical approaches. Konrad [14] first tried to formulate a 3-D FEM 
with three vector components to represent electromagnetic fields in anisotropic but 
loss-free media. The tensors were therefore assumed to be Hermitian in his study. 
A few years later in 1980’s, the number of publications in this subject increased 
typically with applications to waveguide structures. Unfortunately, the variational 
approaches reported by different authors during that period consistently led to non- 
standard and non-Hermitian eigenvalue systems (even with the aid of an adjoint 
system [15,16]). Even worse, the numerical systems derived in this manner were 
usually doubled in size. As indicated in [17]. when a non-standard eigenvalue system 
was manipulated to reduce to the standard form, the size of the system was doubled 



again. Similar r<‘ports \v(m<' al>u seen in laU'r paper." ;l^ lot the pruMeni" ul u.ue 
{propagation inside anisot ropic media. No radiation and "eattering analysi- li< 
rejjorted in this ronte.xt. 


a'- hmi 


In wliat follows, we generalize the FEM formulation to lo.ssy and anisolropu 
electromagnetic problems. Specifically, we show that two different methoils. one wit h 
the aid of an adjoint auxiliary system and the other with the Lagrangi' multiplier, c an 
be used to construct the pertinent functional. Galerkin's method is then compared 
to these two variational techniques. 


2.1.1 Pertinent Functional For Lossy/Anisotropic Media — I 

As is known, the na/ura/ variational functional no longer exists for non-Hermit ian 
operators since no matter what definition is given for inner products (see (2.7) or 
(2.8)). one cannot obtain a self-adjoint operator necessary for vaiural functional 
design. In these cases, we consider a generalized functional 


T =< tF>-<4>,Ka>-<^.K> 


where $ is the unknown solution function of the original PDE problem and K is 
the right-hand-side function as in (2.1). Similarly, is the solution function of the 
adjoint PDE such that 

( 2 . 12 ) 


where can be derived from 

< 'F >=< > 


(2.13) 


with C ^ La- 

It should be remarked that the functional (2.11) reduces to (2.6) (except for a 
possible constant coefficient) if L is self-adjoint. .Also, the original PDE and its 


acljoilil rountrrpart (2.rj| can lx* iliriMicli ilie- \ariai iwn«tl with 

r(‘S|KT! to i1k‘ functions 4> and res]>r('l i\ riy itlu' siiii|)l(‘ d<‘ri\ath)ii uiiiitird 
hore*). After discr(*t izat ion is rarnrd out. tiu' final nunuuical s\>i(‘m i> >\inineiric. 
Tills can be shown as follows. Let 


' j 

where V, is the basis function used for both unknown functions and .i\. i/, are the 
corresponding expansion coefficients. Inserting (2.11) into (2.11) yields 


< CV,.V, > < V,.Ka > - ^ t/, < V,.K > (2.15) 

' ' j 

Ijpon performing the differentiation with respect to x, and yj individually, we get 
the two decoupled systems of linear equations 




;2.16) 


where the matrices and the column vectors K^. K*' are given b\' 


Q- = <£V,.V,> 
Qf, = <£V„V,> 
Kf = < V„ K > 
Ky, = <V„Ka> 


In general, Qy 7 ^ Q", and Q^. However, 

These relations indicate a loss of symmetry of the original problem, but the symmetry 
holds for the overall system! 

The storage requirement is a function of N/2. where .V is the dimension of (2.16). 
Even though there is an auxiliary system needed to complete the analysis, in practice 
this system does not require storage. 



2.1.2 Pertinent Functional For Lossy/Anisotropic Media II 


An alternatixT to using an adjoint syst(Mii is tu cnipluy llic l.aiiranet' inuhi|>lu'i 
tecliiiiquf in constructing the pertimuit functional. The Lagrange inultii.rKM is UMiallv 
used to incorporate additional constraints to a system. To illustrate tlu- teclumpie. 
consider the same PDF model as described in (2.1). and we first rewrite it as 


C<t> - K, = 0 


Next, we assume an expansion space where the solution is defined and solv(>d. The 
unknown function ^ and the multiplier function A are expanded in the same space. 
If (2.17) is regarded as a ‘'constraint", we try to add the constraint to a "null system 

and get 

A) =< A,£$ - K. > 


This functional is now used to formulate the FEM. As described above, on applying 
the Rayleigh-Ritz procedure to both and A using the same set of basis functions. 

viz. 




(2.19) 


we obtain 


A) = Vj,£V,' - K,- > 


( 2 .- 20 ) 


2 3 


Carrying out differentiation with respect to x, and j/j, individually, yields 

i y ^ 


0 Q" 

Qy 0 



( 2 . 21 ) 


where Qf^ = Q], as in (2.16). We observe that (2.21) is similar to (2.16) with two 
decoupled subsystems of the same size. The properties of the subsystems arc also 



similar lo t hosv in ( 2. 1(> l. \\(‘ furl lu'r ol>s('r\ r l hat t 1 k‘ Laui <tm!r muh i t t hnitjiu 

may regard('cl as a special ras<‘ lu (2.1 h wlirif a lioniuunu'uuv adjoint 1M)1. i> 
now virtually assumed (i.c. = 0). .-\oain. it should lx' remiarkt'd that ulu'ii ii 

self-adjoint problem is considered. th(‘ aho\(‘ formulations (both i1k‘ adjoint ‘^\ste^) 
approach and the Lagrange multiplier technicjinv) will result in a symnu'tric s\stem 
and the auxiliary subsystem becomes either rt'dundant (for adjoint approac h) or un 
necessary (for the multiplier technicpie). In the later case, tlu' multiplicu (an Ix' 
considered identical to the unknown function 4> expanded using the same basi> fum - 
lions. This results in an interesting coincidence with Galerkin technique d('scrib(*d 
next. 

2.2 Galerkin Formulation 

Galerkin's method is now considered to formulate the finite element method. Tra- 
ditionally. Galerkin s technique used in conjunction with integral equation employs 
the same testing and expansion functions to obtain a symmetric dense numerical 
system. However, in the case of the FEM, Galerkin's method does not always lead 
to a symmetric system. Apart from boundary conditions, the linear operator of a 
PDE problem determines the symmetry feature of the resulting system. 

Also, unlike the variational approach, Galerkin's method solves the weighted PDE 
by a testing process as 

< >=< ^/.Ki > (2.22) 

w’here $ and ^ are both defined in the same function space. Specifically, in Galerkin's 
method one seeks the solution for the unknowm function ^ w’hich satisfies certain 
prescribed constraints and (2.22), with the aid of another arbitrarily chosen function 



Siniilarlv the varialional a])|)roacli. ilu' Hayl('i'ji,li Kit/ proi rdurc ma\ <il-o hr 
used to projf'ft the continuous si)ace onto a fiinte discrete sej)aral)le lidheit 
Tlie inatliemalical prol)lern is tlieii replirased to seek a discrete solution set uliu'C 
entries are the coefficients of the expansion. Tlie testing function 'I' must, of course, 
be defined in the same discrete space to ensure that tlie original PDL is solved with 
proper boundary conditions. 

It is obvious that if the linear operator £ is self-adjoint, the choice of the testing 
function 'I' = $ results in a symmetric numerical system. Otherwise, no matter how 
the inner product is defined, the final discrete system in general does not exhibit 
symmetry property. 

We observe that Galerkin's method can usually be applied to any linear operators 
even when the corresponding natural functional does not exist. Also in the general 
cases (as considered when describing the functional approach). Galerkin s method 
leads to the same numerical system as the desired portion in (2.16) and (2.21). This 
can be demonstrated as follows. Inserting (2.14) into (2.22), we readily obtain 

E cv, >= E y. < v„ K. > 

i J 


or 


Ek |e^' < 


0 


(2.23) 


.As assumed, ^ is an arbitrary function. Thus the term in the curly bracket should 
vanish, yielding 


Q"x = K’ 


(2.24) 


which is exactly the same as the subsystem derived from the variational approach. 



It is not (“ci that if oiK‘ chose t he Lacratiae mult i|>lier i in the \ ariat luiia! ,i| i| >i nai li i 
as the testing function vj;. tlie same numerical systtun would result. .\ual\ ! u all\, 
though, the entire system is twice ih<* size of that th-rived \ ia Cialerkiii's na'ilu.d. 
This is due to the fundamental clifTereiice of the two teclmi(|ues. He<iai dless. as 
expected, the obtained numerical systems of interest are virtually identical! 

2.3 Total Field and Scattered Field Formulations 

In this section, we focus on a general scattering problem as illustrated in fig. 2.2. 
where a perfectly conducting electric (PEC) body is coated with a dielectric la>ei 
whose relative permeability and permittivity are Jij and fj. res])ectively. (.Note that 



Qrf: dielectric coated region 
fl/: free space = 1) 

Via: absorbing layer (fa,^a) 

Tp: boundary of the PEC body 

Trf: boundary of the dielectric coating and free space region 
T }•. boundary of the absorber and free space region 
To: PEC boundary of the outer absorber 

Figure 2.2: Illustration of a scattering problem setup for scattered field formulation. 

for the purpose of generality, the medium is assumed anisotropic.) 

The situation with absorbing boundary conditions for truncating the FEM do- 
main has been analyzed before (see e.g. [24] or [25]). However, two issues associated 
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uilh tills tv|>e of problems liave not been carefully addresM'd in the liteiatnie. Oiu 
of them is the equivalence between the variational and the (lahukin s method when 
the scattered field is used as the working variable. Proof of this etiuivalnue can be 
tedious and cumbersome but it is nevertheless an important issue. 1 nfortunately. 
one is used to assuming that these two formulations are equivalent without proof. 
.Another issue relates to the recently introduced perfectly matched absorbitig mati'- 
rial. As it turns out. there are several advantages to use artificial absorbing materials, 
including accuracy control, conformality, ease of boundary treatment, etc. Ho\\e\ei. 
their inclusion introduces additional artificial conditions inside the absorber layer 
and care must be taken when those conditions are enforced in the FEM formulation. 
Moreover, although a metal-backed absorber layer simplifies the FEM implementa- 
tion, the multilayered FEM region contains high inhomogeneity, which again requires 
a careful presentation of the formulation. To this end, we extend our theoretical dis- 
cussions on the FEM to scattered field representations, where the treatment of the 
boundary and transition conditions will also be described. 

2.3.1 Scattered/Incident Fields and Boundary Conditions 

Referring to fig. 2.2, we begin with the wave equation in terms of H (the E 
formulation may be readily handled by duality). To proceed with Galerkin's method, 
we first write as 


pjtot jjscai 


(2.25) 


where and H'"" are the scattered and incident field, respectively. Next, weight- 

ing the source free wave equation with the testing function V yields 

/ V ■ { V X !;' • V X H - klfr • h} dQ = 0 (2.26) 

Jq 


\vhf*re n — Hi + ilf Ha eiicom[)ass(‘s tiio ciitirr ctMiiiMitat u)iial tviiion. lu |uuc I'nl 
with tlif“ cieri\ation of the weak form wa\'e <“(|uaiK)ii. it is ner(‘ss<ti \ tii iiii tudurr 


certain constraints on the scattered and incident fields within H and on the l)onndar\. 
First, since the incident field is iiol allowed to pass through the ahsorhei la\cr as 
well as the metal back wall Fq. we note that 


H'"'(r) = 


pj5Caf ^ 


H 


scat 


r^Od + Qj 

r € otherwise 


f) -K 


with the incident field satisfying 
{V X • V X 

It is thus evident that the scattered field satisfies the homogeneous wave equation in 
regions Hq + Qj and the inhomogeneous wave equation in flrf. 

The boundary conditions on can be readily derived by consistentlv applying 
the field decomposition (2.25). Note that in accordance with (2.27), the electric field 
is likewise decomposed as 


= 0 
^0 


r € rio + r?/ 
r G otherwise 


(2.28) 


= E*"“' -f E’"" (2.29) 

However, one should be cautioned that E'"'^ and do not satisfy Maxwell equa- 
tions in the dielectric region. That is, 

ju;eoE'"‘^ ^ y X reQd (2..30) 

which conflicts with what one would intuitively assume. Conventionally, the incident 
field is assumed to exist in the dielectrically coated region Clj as if there was no 
dielectric there. After a quick glance, one would immediately arrive at a conclusion 
that (2.30) is against Maxwell theory. In reality, it can be proved that if E”‘" and 





H”‘' iii(k‘<‘d satisfied Maxwell etjuations in n,;. a contradicit)i v l>u»mdarv ( oiuiitioii 
would immediately result. This can he seen 1)V taking a curl operation t>n holh suh" 
of (2.25). In view of (2.29) and .Maxwell equations in dielectric- media. w«' would have 


r' • V X H'"' =r‘ • V X H 


— 1 


■ V X H*-'"' r € ilj 


( 2 . 51 ) 


Imposing the condition of total field tangential continuity at the boundary I ,j woulc 
yield 


X {?;■ . V X H-'-'lr; - f,-‘ ■ V X = 0 (-'.321 

which is a homogeneous Neumann boundary condition. However, if we start with 
(2.27), upon taking the curl operation and imposing the condition of total field 
tangential continuity at the boundary Tj, we get 


n X 




. V X H 


scat 1 




zz-l 

- V 


• V X = -n X {c/ -f/} • V 


X H' 


( 2 . 33 ) 


which is an inhomogeneous Neumann boundary condition. 

This inconsistency is because (2.32) was derived on the basis of the decomposi- 
tion (2.29) and the assumption that the incident field within dielectric regions also 
satisfies Maxwell equations. However, the decomposition (2.29) is artificial and it 
is therefore necessary to keep in mind that only (2.27) holds true when deriving 
boundary conditions. 

As a rule of thumb, an appropriate interpretation of the phenomenon should read; 
the incident field inside dielectric media existed in the same fashion as in free space 
as if the free space was replaced with the media. Mathematically, this implies the 

condition 




;£oE'"^ = ly • V X H 


r 


(2.34) 



:n 

( onsisteiitly one can derive tlx* otlier lioundary romiitioii> re(|iiirr<l fui ilir I I'M 
formulation following the same proeedure. They are classified as Dirichlel am! Neu- 
mann conditions as follows. 


• Dirichlet Conditions 


Boundaries 

Conditions 

Tp 

7? X = K° — 7? X (Ks unknown current) 

r, 

n X = 0 

T; 

n X = -75 X 

To 

h X = Kj (Kj unknown current) 


• Neumann Conditions 


Boundaries 

Conditions 

Tp 

n xz/ - V X H*"“' = - 

“77 X ' V >< 

H>nc 



-1- 

+ 


n X |f/ • V X 

=-l 

= — n X 

•V X H'"" 



+ 


T/ 

n X |e/ • V X f 

~i - 

= n X ^ 

C X H'"" 

To 

n X e/ ■ V X = 0 



where {} |^ denotes and and should take the values at 

positive and negative sides of a specific boundarv. 

2.3.2 Galerkin’s Method 

Returning to (2.26) and in view of the vector identity 

A • V X B = (V X A) • B - V • (A X B) (2.35) 

and the divergence theorem, we obtain the corresponding weak form wave equation 


Intd + Inty + Into = 0 


(2.36) 






where 




V X V • V ,> (H-'-''-' 


H"‘ ) - A--;V • (H' 


(/n 


+ 


Int 


J V • X • V X h) (/>■ + V • X 7,;' • V • H”' ) </> 

7 / V - (n, x7;’ ■ V X ^/5 

= f V X V ■ fy' ■ V X (in-\- f V ■ ^-ni X fj' - V X H'"") (/< 

Ju, 

+ / V ■ (h2 X !;' ■ V X </5' 

Int, = f V X V • !;' • V X H*^“' (/n + / V ■ (-»2 X f;' • v x ds 

Jua 

+ / V • (ns X !:' • V X c/5 (2.39) 


with no, ni- ^2 and ns being the unit normals at the boundaries Fp, Tj, T; and Tq. 
respectively. They all point away from the center of the PEC body (i.e. outwards). 
Invoking the boundary conditions as tabulated above, we observe that the surface 
integrals on Ep in (2.37) and on To in (2.39) vanish. Also, the sum of the surface 
integrals on P; in (2.38) and (2.39) reduces to 

- f V-n2xf7'-VxH*""d5 (2.40) 

Jr, 

Similarly, the sum of the surface integrals (involving on Ej in (2.37) and (2.38) 

becomes 


I V-n.x (f;‘-f;‘) vxH' 


dS 


(2.41) 


Note that both integrals (2.40) and (2.41) will not contribute to the system, but to 
the excitation on the right hand side. Of course, the excitation term (2.41) would 
disappear if the permitivity is continuous across the boundary E^. By duality, for 


an electric field formulation, a term similar to (2.41) will appear for discontinuous 
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l'hi> observation do<“s not liuld for CJ.tOi. uhieli ari-'C' fiom tin' inlu-i 
ent incident field definition on the two sides of the lH)iindar\ I a> niai lieinat u al ly 
shown in ('2.27 ). Furthermore, tlie second term in (2.11 ) will l«- (ancelled uiu with 
the integral on F^ involving H'"" in (2.37). .After a simi)le manipulation, the tinal 
svstem reduces to 


V X V ■ t, ' • V X - k^V ■ T^r ■ (19. 

= - / (v X V • !;* • V X H”"' - kl\ • ■ H”"') (19 

-f V -n, X !;' • V X H'""d5+ / V-n 2 

Jr, 


X (/ . V X H""' ds- (2.42) 


where Cr is again the relative permittivity with respect to the specific region. This 
is the desired weak form of the wave equation, and a similar equation was obtained 
in [24] using the functional formulation for isotropic media. 


2.3.3 Variational Method 


It is now of interest to employ the functional formulation to obtain the equation 
corresponding to (2.42). To this end, it is intuitive to begin with the total field 
representation and use the functional 

= f (v X H-f;^ ■ V X -h) df) (2.43) 

“ Jn^+n,+Qa ^ ’ 

where, as before, Ir and denote the corresponding relative permittivity and per- 

mebility, respectively. Also similar to Galerkin’s method, one would now logically 

proceed with the field decompostion H = and express the functional 

in terms of the scattered field. It is unfortunate that this approach wdll result in a 

different form of linear system than (2.42) and the detailed mathematical proof is 

presented in Appendix D. 



1 lie <.lisrr<‘[)ft iK"\ hi' isos from t he hssuiii|>t u)ii ol l li<' fuiict loiihl i 1? i . u li ii 1 1 i' i u >i 
a valid (’Xi>rc>sion wlion tlio corresi>ondiiig I’DK oixTator is no 1oiic<t srlf-mlioini Ko- 
call in section 2.1.1 that the preseiic<‘ of general anisotroi)i( and lossy ineilia reiiuiies 
the application of an au.xilary adjoint system, together with which the pi'rtmeni 
functional is given b\ a generalized form in (2.11). It is therefore necessary to begin 
with this generalized functional rather than (2.43). In a source-lree region. (2.11) 
explicitly takes the form 

.T(H,. H) = H, • { V X • V X H - • h} r/\ (2.44 ) 

and the variation is imposed on the adjoint variable Hq to get 


6.T(Ha,H) 


= 0 

«H=o 


(2.45) 


The variational functional method of this version proves valid and identical to Galerkin 
technique for any linear operators in electromagnetic problems. It is also interesting 
to note that once compared to Galerkin’s method, the adjoint field quantity Ha in 
(2.44) seems to take place of the testing function V in (2.11). However. Ha theo- 
retically differs from V in that the former is defined as a solution function for the 
adjoint system of the original PDE problem, whereas the latter is just an arbitrary 
admissible testing function that does not have to be a solution to the adjoint system. 

.\part from the concept difference as mentioned above, the mathematical proce- 
dure required to derive the FEM system is similar to that presented in the previous 
subsection, and will not be discussed here to avoid repetition. One can be assured 
that the final system obtained from (2.44) and (2.45) is of course identical to (2.42). 



2.4 Parameter Extraction 


An acrurato full wavr analysi' can only prc'dict the lu-ar licM i f'ui l’l)K t> i'f 
methods) or current (for integral based teclmi(lue^i distributions, which can be used 
to obtain certain practical parameters dej)ending on applications. For intricate sys- 
tems. involved numerical models may be needed for outinit data e.xtraction. iiic lucliug 
far field evaluation in the presence of non-canonical platforms and a de-embedding 
process for antenna feed network or circuit simulations. 

.Antenna parameters can be readily evaluated after the near field distribution 
is achieved via full wave analysis. The de-embedding process is required for feed 
network or circuit modeling and will be discussed in chapter 7. For a non-planar 
platform the far field evaluation can be obtained from the general discussion in 
chapter 1, where the formulation in terms of the dyadic Green's function must be 
used to consider the equivalent currents and the free space Green's function. 

2.4.1 Radiation and RCS Pattern 

In the case of antennas, we are mostly interested in their radiation and scattering 
patterns and other related parameters such as gain and axial ratio. (The near field 
quantities such as input impedance, feedline S-parameters, etc., will be disccussed 
in later chapters). Both radiation and radar cross section (RCS) patterns can be 
readily characterized with respect to the 3-D spherical coordinates 6 and 0. 

Consider the planar cavity-backed antenna as showm in fig. 2.1. Once the field 
distribution on the aperture S of the conformal antenna is obtained from the full wave 
analysis, we can then proceed to evaluate the far field pattern. The most straightfor- 
ward approach for this computation in the presence of an infinite conducting ground 
plane is use of the equivalence principle. To do so. we define the magnetic current 





as M = '< E'. where E'’ is the elertrie field (waluatetl on the apevnne ami : i' 

normal to the aperture surface. Tli<‘ eh'Ctrie vector potential i- then given 1>\ 

-jkr 


F(d.O) 




C-’.-K' 


where the electric field is typically expanded in terms of the surface vectiu l>asi^ 
function S'. Introducing this expansion. (2.46) becomes 


F{0,o) = 


C{)€ 


-jkr 


4nr 


V 


( 2 . 47 ) 


where V' denotes the sum of the surface integral over each of the discretization 
elements on the aperture. The far zone magnetic field H becomes 


Hj -- —jutFr 


( 2 . 48 ) 


where Hr represents the transverse component of the far zone magnetic field, whose 
0 and 4> projections are given by 


He = -3 


. toe 


-jkr 
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-jkr 


(2.49) 
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vhich Pe = -j{-] and P^ = and {•} stands for the corresponding terms 

in the curly brackets. The RCS of the t {t = 0 or <p) component can be represented 
in terms of P$ and P^, to yield 
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(2.50) 



1(1 


where A and Zo = are llie fre<‘ space wa\<'lemitli and t he ini rin>ic imped. incr. 

respectively. In '>0|. the incident wave|//'| was s<‘t to unity witlioul loss ol uem i al 
ity. In })ractice. it is customary to ('xi)ress lU'S in <1B. where rr. is usuall\ nonii.di/ed 
to A* or to squared meter first. 

The radiation pattern may also be represented in the same manner as mentioiu'd 
above. The difference in procedure here is tlie normalization with respect to a max- 
imum radiation field value. The reason is that in radiation mode, the antenna is 
excited by an interior source rather than the incident plane wave //' as in tin- scat- 
tering case. Of interest therefore is the relative field intensity in the far zone. To get 
this, w'e represent the far field intensity in terms of the above calculated quantity 
cr’’'® to avoid a repetition of post processing. Specifically, the formula 


= 




(2.51 


is used for radiation analysis. 

2.4.2 Gain and Axial Ratio 

Gain (G) and Axial Ratio (AR) are two important parameters which indicate the 
antenna’s performance. It is also noted that these two quantities typically charac- 
terize the far zone features of the antenna. 

By definition, the gain G for a lossless antenna (w’ith 100% efficiency) is given by 




2ttU„ 


rad 


(2.52) 


for a cavity-backed structure, where T'max is the maximum power density and Prad 
denotes the total radiated power from the antenna in the upper half space. (It is 
noted that this definition of G is identical to that of directivity for lossless antennas.) 
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Since tlie maNiimim power dcMisity I can al>o l>t‘ expresx’d a- 


( „,a^ = — 
4 <i 


( J. >■> ' 


tlie antenna gain heconicb 


G = 


Zo{(Te + f^c 

2Pr.d 


Thus, the computation of G is rather trivially done once a is found as given in the 
previous subsection. In reality. P^ad may be evaluated on an assumiition that all 
input at the feed is transferred to antenna element(s) and radiated. In this case, one 
has Prad = ^^Pin- ''here / is the known current source on the feed, and P,,, is the 
input resistance of the antenna measured at the reference plane. 

However, this scheme may not always work since the gain (more precisely dir(c- 
tivity) reflects the far field behavior, while the input resistance computation relies 
on an accurate full wave model for near field prediction. It is well known that the 
far field and the near field computations offer different accuracy. This accuracy in- 
consistency arises if Pin is used to determine the gain G, a far zone pattern whose 
accuracy is directly governed by the near field computation without averaging effect. 
To avoid this accuracy inconsistency, one may calculate the gain or the directivity by 
evaluating Prad from the far field radiation pattern. Specifically. Prad can be obtained 
by integrating the radiation intensity 


Prad 



( 2 . 55 ) 


over the half space. It is obvious that if a certain symmetry of the pattern remains 
such as circular about the vertical Z axis, then t'(6>,d)) reduces to P(e) and the above 


IJ 


integration can be efb'ctively evahiat(‘d. Ot lierwi^e. a nuiiieric.il 1) iiiteei <it urn i- 
required. 

-•\xial ratio (,4/f’) is another important antenna paraim'tt'r. esp(>ciitll\ when a 
circular polarization (CP) of antenna s performance is of the |)iimar\ concern. .'sinc<‘ 
.AR also features the far zone effect of the antenna, it is desirabU' to determiiK' Al\ 
with a minimum computational load. It is noted that given the above jire-calculated 
ae and cr^. one is again able to determine AR uniquely by 


AR=^ 

' short 




with 


Viong = Y 9 { + CT 0 + 2 o -|(72 cos 0 ] " I 


Vshort = y 9 - [(^9 + 2a|cr2 cos /5] " | 


(2.oT) 


where 0 = 2($//^ — the twice of the phase difference between the two magnetic 

field components. It can be obtained from the quantities Pe and defined in 
(2.50). Since these two quantities are complex numbers, the phase difference is 
readily represented as 


0 = -j Im 



(2.58) 


with Im {•} being the imaginary part. 



CHAPTER III 


Edge-Based FE-BI Technique 


3.1 Introduction 

Numerical methods have been serving the engineers and researchers for many 
years in antenna analysis and design. Among them, the moment method in con- 
junction with various integral equation (IE) formulations played a major role [1 I]- 
However, IE methods are associated with field representations in which the appro- 
priate Green's function for the specific geometry must be employed and this limits 
their versatility. Moreover, IE techniques are usually formulated on the assumption 
of an infinite layered (not inhomogeneous) substrate, a model which deviates from 
the practical configuration and leads to inaccuracies for larger bandwidth antennas. 
Furthermore, in the context of IE methods, antenna excitations are represented us- 
ing simplified models that differ more or less from the actual configurations. .Also, 
due to the singularity of the current distribution near the feed junction(s). special 
measures must be taken [26] for proper modeling. In contrast, the hybrid Finite 
Element-Boundary Integral (FE-BI) technique alleviates these difficulties and this 
was recently demonstrated when the method was applied to inhomogeneous objects 
of canonical shape scattering [27,28] and rectangular patch antennas [9]. 

Based on the past success of FE-BI methods for antenna analysis, it is desirable 
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to cxIcikI t 1i«- nu‘t IkkI to antennas of arbit rar\ slia|x'. In t lii^ clia j>t ei . an talcc b.i-ed 
liyl)rid finite eieinent -t>oiindar\- integral fonmilation is jireseiited fui the (liaiaeiei- 
izatiori of arbitrarily shaped cavity-back('d antennas [L’h;. An example of vui li a 
configuration is shown in fig. .‘3.1. where a cavity is recessetl in a metallit. tuoumi 
plane enclosing the FEM volume. The antenna elements on th<‘ aperture ma\ be 
excited by different schemes, such as a simple probe, a magnetic frill generator, a 
practical coaxial cable, microstrip line, slot or a CPW line. In the context of tht' 
FEM. the cavity is first discretized into a number of tetrahedral elements that natu- 
rally reduce to triangles on the cavity's aperture. For non-rectangular patches this 
triangular gridding is. in general, non-uniform and the exact boundary integral for- 
mulation based upon this mesh applies to any patch shape. As a result, the hybrid 
FE-BI technique is capable of modeling arbitrarily shaped cavity-backed antenna 
configurations, different substrate inhomogeneities, anisotropies, as well as various 
practical excitation schemes. 

3.2 Hybrid System Functional 

In this section, the edge-based hybrid FE-BI method will be formulated using 
the variational principle, where matrix algebra notation is employed so that one can 
readily extend the formulae to the general anisotropic case. As presented in [9]. 
the complete functional pertinent to the scattering and radiation by a cavity-backed 
configuration shown in fig. 3.1 may be written as 
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Figure 3.1: Illustration of a typical radiation and scattering problem. 


F(E) = |(V X E)- X E) -fcotrE-Ej c?r 

+ 2jkoZo jj [E xH‘)- ~dS 
+ E • ^jAroZoJ. + V X dv 

-2kl yy(Exi)-|yy (Exf)-(^I+^W^G'o(r.r')<i5'|c/5 


(3.1) 


where J, and M. represent interior electric and magnetic current sources within 
the cavity V; H' is the incident field, if any, from the exterior region; the surface 
S encompasses the cavity aperture excluding the portion occupied by the antenna 
elements; Cr and pr denote, respectively, the relative permittivity and permeability; 
ko is the free space wave number, I the unit dyad, and Go(r, r') the free space Green’s 
function with r and r' denoting the observation and integration points. 


3.2.1 FEM Subsystem 


In proceeding with the discretization of (3.1), it is convenient to decompose it as 


F = Fv + Fs 


(3.2) 



I(. 


where f \ denotes the \ ulume integral cunt rihnt iuii' ami similarly / . actuimt- for t lie 
surface intt'gral contributions. The cavity \olume i> sub(li\ idc-d into N t et r.iliedi al 
elements (e=l.d — ,\). and within each tetrahedron as shown in lig. d.'J the iield 

1 



Figure 3.2. A tetrahedron and its local node/edge numbering scheme 


is expanded using edge- based elements as 


with 


e = [v]J{e}, 


[V'le = [{v^x}{Vv}{v;}], 




{v;} = 


^u 2 


u = x,y,. 


\ 
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Vu6 
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[E]. = 


E 2 


\EeJ 


(3.4) 


e 



in whirli \ is the a {u = .t.ij or r) coini)uiK'iil of the volume vt>ctoi l>asi- iumiiuii' 
along the /th edge. The unknown sector {A }, ha> six entri('s. i>ue for em h letrahe 
dron edge. (Here we use square brackets for matrices ami curls biacki't'- hn scilomi. 
Inserting (3.3) into (3.1). and taking the first variation of F\ with resiiect to {/.},. 
yields 

^Fv = + {A'}.} 

where 


[.4]e = JJj^ !^^[DVUDV]l (3.G) 



To carry out the above integrations, it remains to introduce the volume expansion 
or shape functions Ve- For our implementation we employed the linear edge- based 
shape functions for tetrahedral elements given in [30,31]. The explicit finite element 
matrix entries associated with a typical tetrahedron (as shown in fig. 3.2) are given 
in .Appendix A for reference. 
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3.2.2 Boundary Integral Subsystem 

To discretize the surface integrals in ( .'M ). tin- a|)(Ml ure i> >ul>ili \ iiled iut^ 1 1 i.iii 
gular elements since these correspond to the faces of the tetraliedrals. Within eac li 
triangle, the field is represented as 


e = [S]J{e,}, 


(d.h) 


where 


:s'l. = ({sn(s.)|. 


/, \ 




^u2 


V / 


{£,}e = 


-s2 


Es 3 


u = x.y 


(3.10) 


\ / e 

in which S^, is the u{u = x,y) component of the surface vector basis functions along 
the ith edge. On substituting (3.9) into the surface integrals of (3.1) and taking the 
first variation of Fs with respect to {£’s}^, we obtain 


= ^ {[»].{£.}. + {i }4 


where 

[B]e = 



-2kl[Se][S,]^ + 2 

Go(r,r') dSdS' 




(3.11) 




dx 


and 


{L],=j2k^Z^jj^\S.] 


HI 


-HI 


dS 
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Note that in (3.12) the elcnu'nts of tin- array i,<! ar<' function-- of the ol.M-rvaiiun 
vector r. whereas the elements of ar<- with respect to the inle|iraliun pomi r'. 


A suitalrle set of linear edge-based surface basis functit)ns is 


S,(r) = 


X (r - r, )<(r) r t .S 


(3.1-n 


0 otherwise 

In this expression (referring to fig. 3.3). /, denotes the lengtli of the /lh edge and r, is 

\fh edge 



Figure 3.3: Pair of triangles sharing the fth edge 


the position vector of the vertex opposite to the fth edge. Since each edge shares two 
triangles, one is defined as the plus and the other as the minus triangle. Therefore. 
c(r) is given by 


where S, = T S; . The constant Ae in (3.14) denotes the area of the plus or 
minus triangle depending on whether r G 5+ or r G 5^. We note that S,(r) x c 
yields the basis functions used by Rao, et al. [32] in their moment method solution 
of boundary integral equations. The explicit expressions for the boundary integral 
equation subsystem is given in Appendix B for reference. 


3.2.3 Combined FE-BI System 


To construct tljc final syslen] for th(‘ solution of tlx* (‘hnUric fu'ld ^ \w 

combine (3.5) and (3.11 ). and after assembl\ \v(‘ obtain tin* svsteni 

{[A]{E} + {A'}} + + {/,}} = 0 (.IK.l 

In this. {A} and {1} are the excitation vectors due to the interior current soukcs 
and the exterior excitation, respectively. The unknown electric held vector {A} 
consists of all held expansion coefficients with respect to the element edges except 
those coinciding with perfectly electrically conducting (PEC) walls. PEC antenna 
element(s) or PEC pins inside the cavity. Finally, the vector {E^} represents the 
unknown surface helds whose entries are part of those in {E) with their corresponding 
edges on the aperture. The explicit expressions for the matrices and vectors in (3.16) 
can be readily extracted from (3.6), (3.7) and (3.12) (see also [.33]). It is evident that 
[A] and [B] are symmetric as a result of the assumed isotropic medium and reciprocit}'. 
In addition, [A] exhibits high sparsity due to the FEM formulation whereas [B] is 
fully populated. 

Two approaches may be followed in carrying out the solution of the combined 
subsystems when an iterative solver is employed such a.s the biconjugate gradient 
(BiCG) method [34]. These approaches differ in the manner used for the evaluation of 
matrix- vector products called for in the iteration steps. One could sum the coefficient 
matrices [.A] and [B] by adding up the corresponding matrix entries prior to the 
execution of the BiCG algorithm, or instead the resulting vectors may be summed 
after carrying out the individual matrix-vector products. We observe that the first 
approach is more efficient in terms of computation time after reordering the combined 
matrix and storing only the non-zero elements. This is because, in the context of 



this sclK‘m<*. tlir coinbinalioii of llu* two matiic'c's is ju'rionnrd uiil\ oikt (.nit>idr \ \\c 
iteration. 11 owover. the second aijproach cuiupat iMe wii li the Bi( (i-l 1 1 ■'i lieiiu'. 
where the 1 FT algoritliiu is employed to ('xploit the coiivoliilional |)ro|)eit\ of the 
integral operator, thus eliminating a need to explicitly store the entire HI matrix. 
The BiCG-FFT technicjue will be discussed in chapter -4. 

3.3 Numerical Implementation 

Based on the above presented FE-Bl formulation, the hybrid method was im- 
plemented and a computer program was developed for the analysis of radiation and 
scattering by cavity-backed patch antennas of arbitrary shape. The antenna geom- 
etry/mesh is first generated as shown in fig. 3.4 and supplied to this program in an 



Figure 3.4: A typical geometry/mesh for a cavity-backed circular patch antenna. 

input file which, as a minimum, contains 

(a) the nodes and their {x,y,z) coordinates; 

(b) the tetrahedral elements and the corresponding nodes forming each element: 


(c) the nodes identifying the cavity aperture; 


( (1 1 I he nodes ident ihing iiH'tallir l>oiimlari«'s. ainemia palclncM. fccdi > i. and i>u' 
sil>lc \<‘rt iral posts. 

For arbitrary antenna geometries, it is necessary to employ a sophist ii at cd Milimic 
mesh generation package and a number of tliese are availal>le commerciall\ . I x picalls 
eacli of these packages generates a "universal file" that can be readilx i)re])i o( ess('d to 
e.xtract the aforementioned input list. Given the above list of data, an interpretation 
routine is used to convert the information from node-elements to edge-elements. We 
usually refer this procedure as data preprocessing. The flow chart shown in fig. 
describes the major implementation procedures from the mesh generation, a few data 
preprocessors, the FE-BI kernel, to the BiCG solution and finally the data output. 

One of the primary complications in the hybrid technique implementation is tlie 
efficient treatment of combining the two separate subsystems. It is noted that the 
FEM sparse matrix is large in dimension but requires less storage space, while the 
boundary integral subsystem is always small in size but can be dominant in terms 
of memory demand. This is particularly true when the non-metallic portion on the 
cavity aperture predominates over the antenna radiating elements. Furthermore, 
the boundary integral subsystem in a general purpose hybrid FE-BI implementation 
is entirely independent from that of the FEM and even the basis functions can 
be independently developed. This also accounts for the two arbitrary numbering 
systems and combining them is a relatively complex task. One major advantage of 
the method however is that these two subsystems can be developed and validated 
individually. 

Once both of the subsystems are verified, the coupling of the subsystems is ac- 
complished by enforcing the boundary conditions implicitly on tangential H fields via 
the integral representation and explicitly on tangential E fields over the interior and 





IMPLEMENTATION FLOW CHART 




Figure 3.5: A flow chart describes the major implementation procedures from the 
mesh generation, a few data preprocessors, the FE-Bl kernel, to the 
BiCG solution and finally the data output. 
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f'Xtcnor r(‘gi<jiis, A signifiCfiiit ('fTort \vii> cl('\\)i<‘d to (l<’\<'k)|)nig ti ])ioiii<im m ■'Ui h <i 
inaiiiirr so tliai f)Otli tho storage and computational re(|uir<-ment '• can la miiiinii/cil. 
Specifically, the boundary conditions on tlie metallic siirfacf- are eidoncd n pnon 
to condense the system which involves only nonzero field components. To this point, 
the sparse finite element matrix was stored as a single array of length .V, wh<-U’ 
.Nt is the total number of unknowns within the ca\it\ \olume and rlenotes the 
maximum number of nonzero row entries. The B1 matrix was stored in diflerc'iit 
\\a\s foi the e\aluation of the matrix-vector products. If the BiC'Cl solution was to 
be carried out without special treatments (such as incorporating the FFT). then the 
As X As BI integral matrix is added to the FE array resulting in a 1-D array about 
^fAnj + A'j^ long. For slot antennas, including cavity-backed spirals, and moderately 
sized systems, it was found preferable to use this scheme. In that case the generation 
of a single combined FE-BI matrix before execution of the BiCG algorithm reduces 
the computational requirements. This is because a number of operations associated 
with the repeated combinations of the FE and BI subsystems within the BiCG iter- 
ation is avoided. The alternative is to carry out the matrix-vector products for the 
FEM and BI subsystems separately. The advantage of the scheme becomes appar- 
ent when a special treatment is performed on to the numerical system for efficiency 
consideration and this wdll be investigated at certain depth in chapter 5. 

3.4 Selected Numerical Results 

We present below^ some representative numerical results for the purpose of vali- 
dating and demonstrating the robustness of the tetrahedral formulation for scattering 
and radiation by different configurations of cavity-backed antennas. In each case the 
computed results via the FE-BI method are compared with reference measured or 
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calculated data. 

Scattering and radiation by a circular patch: 

Fig. 3.(. illustrates a circular patch rcsicling on the surface of a U. ID(- cm thick 
substrate having a relative dielectric constant of = -M'- 1 he patch s diameter is 

2.6 cm and the substrate is enclosed in a circular cavity 6.292 cm wide. 1 his cavity 
and the patch were recessed in a low cross section body for measuring its H( S. .\ 
comparison of the measured and calculated backscatter ago H(’S as a function of 
frequency is also shown in fig. 3.6. For this computation the direction of the iiicich'iit 
plane wave was 60° from normal, and as seen the agreement between measurements 
and calculations is very good throughout the 4-9 GHz band. Input impedance mea- 
surements and calculations for the same patch are displayed in fig. 3.i. The probe 
feed in this case was placed 0.8 cm from the patch's center and it is again seen that 
the measurements and calculations are in good agreement. 


Radiation by a one-arm conical spiral: 

We considered the modeling of this radiator to demonstrate the geometrical ver- 
satility of the FE-BI method. Two projections of the spiral radiator and surface mesh 
are illustrated in fig. 3.8. The top and bottom edges of the strip forming the spiral 
follow the lines p = O.O5O3Aexp[O.221(0±2.66)], r = a± exp(0.221d>). where (p.O.c) 
denote the standard cylindrical coordinates, a± are equal to 0.0832A and 0.0257A. 
respectively, and 0 < <^ < 27 t. This spiral arm resides on an inverted cone (9.24 cm 
tall) whose bottom cross section has a diameter of 1.68 cm and the top cross section 
has a diameter of 21.78 cm. For our calculations A = 30 cm (/ = 1 GHz) and the 
spiral was situated in a circular cavity 10.01 cm deep. The computed principal 
plane radiation pattern taken in the <f> = 90°-plane. using a probe feed at the ca\ it\ 
base, is given in fig. 3.9. It is seen that this pattern is in good agreement with the 
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Figure 3.6: Comparison of the computed and measured aee backscattei RCS as a 
function of frequency for the shown circular patch. The incidence angh' 
was 30° off the ground plane. 
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Figure 3.7: Comparison of the computed and measured input impedance for the 
circular patch shown in fig. 3.6. The feed was placed 0.8 cm from the 
center of the patch and the frequency was swept from 3 to 3.8 GHz. 
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data ojvcn in As can l)c expected, the /A pattern iiiui >lu>\vni ilitfered Ituni 

the measured data near the horiiton i>ecause of interh’ience from the linite ( liculai 
cavitv housing the spiral wliicli was included in the analytical model. I he laltei ua> 
not part of the measurement configuration which consisted of the spiral antenna on 


a large circular plate. 

Annular slot impedance: 

Fig. 3.10 shows a narrow circular (0.75 cm wide) annular slot situated m a cir- 
cular cavity 24.7 cm wide and 3 cm deep. Because the annular slot is narrow, the 
implementation of the BI subsystem is very small for this application and as a result 
there is no need to invoke the FFT in the BiCG algorithm. The FE-Bl method is 
basically quite effective in modeling small aperture configurations without a need for 
special computational considerations. Input impedance calculations as a function of 
frequencv for this radiator, excited by a probe placed across the slot, are shown in 
fig. 3.10. and agree well with the values calculated via a modal-boundary integral 
method [14]. For these calculations, the frequency was swept from 700-1000 MHz. 
The dielectric constant of the material filling the cavity was set to tr = 1 -35 as in [36] 
and this is an effective value to account for the presence of a dielectric slot cover used 
as part of the measurement model for holding the plate. 

Stacked circular patch antenna: 

To demonstrate the capability of the developed hybrid technique, we now present 
a qualitative study and visualization of the near field distribution inside a cavity- 
backed, stacked circular patch antenna as shown in fig. 3.11. Note that the similar 
configuration with rectangular patch shape has been investigated and found a sig- 
nificant bandw'idth increase. This is because of the dual frequency resonance due 
to the tw^o stacked patches. The circular patches are more attractive than stacked 
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Figure 3.8: Illustration of the configuration and mesh of the one-arm conical spiral 
used for the computation of fig. 3.9. 



Figure 3.9: Comparison of the calculated radiation pattern (£’^), taken in the o = 
90°-plane, with data in reference [35] for the one-arm conical spiral shown 
in fig. 3.8 





d = 3 cm 





Figure 3. lOr Comparison of input impedance calculations for the illustrated cavit 
backed slot. 
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rcctaiiffular patrlu's Ikthusc they orrupv a Miiall area u licii opcralcil at l lie >amc lie 
(jiu’iiry [-5(]. 1 nfort unatel\ . no .siifiiricnt rosoau li on llii> ‘irunu'lrv lia> Ixvn ivpui iril 
in tlir literature due to a lack of analysis tool. 

It turns out that the above presented hybrid FEM techni(iu<' is well suitt'd for 
this study. To show this, we chose a cavity-backed, stacked circular j>atch antc'nna 
fed with an offset single vertical post underneath the lower j)atch to link the ca\it\ 
base (viewed as a ground plane). No direct electric contact between the u|>per and 
the lower patches exists and thus the power transfer has to rely conij)letely on the 
electromagnetic coupling from the lower to the upper patch. This can be clearlv 
verified from the near field visualization, which is available and complete only via the 
PDE related techniques. (The laboratory measurement may provide the image of the 
field distribution above and near a microstrip [3S].) Another interesting point is that 
though the antenna was fed with a single offset probe, the energy is concentrated at 
two distinct regions. One is around the probe feed, and the other is near the opposite 
location with respect to the center of the patch. The two regions act as out-of-phase 
electric pole to effectively excite the antenna. Although the patches are circular in 
shape, the offset excitation ensures the linear polarization in radiated fields. 



Figure 3.11: Visualization of the near field distribution at the lower layer of a stacked 
circular patch antenna. 



CHAPTER IV 


Efficient Boundary Integral Subsystem — I 

4.1 Introduction 

As is known, the hybrid finite element-boundary integral method is accurate and 
capable of handling a variety of conformal antennas. However, the drawback as- 
sociated wdth this technique and any other global truncation approaches can make 
it less attractive. This is especially true if one is interested in modeling large an- 
tenna svstems (arrays). Although the FE-BI method is particularly suited for the 
configurations with relatively small aperture size and possibly complex cavity design 
(including feedlines, isotropy/anisotropy, other layers of metals, etc.), it would be 
much more useful to accelerate the speed and reduce the CPU requirement for the 
hybrid approach. One possible solution is the CG-FFT technique discussed in this 
chapter. 

The boundary integral (BI) equation subsystem leads to a fully populated matrix 
whose size is determined by the number of aperture mesh edges. For large apertures, 
this analysis becomes impractical in terms of storage and computation time require- 
ments, and to overcome this inefficiency, a uniform zoning of the aperture is required. 
By resorting to the structured mesh, the boundary integral matrix can be cast into 
a discrete convolutional form, thus permitting the computation of the matrix- vectoi 
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|>r(Kliicts \ ia the discrrtr Four'u'r traiisforni (Di l i and au)idiii<i a invd lu 'luic llir 
full B1 matrix. This memory savimi scIk-iih* ha.' alr<'ad\ Ixvii a))))lit‘d tu 11, -uliinuii' 
involving rectangular grids and a similar implenuMilation was al'ii ic|)oited 

for triangular surface grids [-10] involving inherent a[)()io.\imat ions. In this cha|)l('i. 
we first show that the BiCG-FFT solver can be precisely implemented on unifoini 
triangular meshes. The differences between the rectangular and triangular meslies 
are also described. For non-rectangular antenna geometries, a special treatment k'- 
ferred to as the overlaying scheme is proposed and discussed in section 4.3. .-\ few 
results are presented which demonstrate the method's validity. 

4.2 Application of Conjugate Gradient Algorithms 

The Conjugate Gradient (CG) iterative solution of linear systems of equations has 
been extensively investigated and the representative references are collected in [41]. 
Although the state-of-the-art CG algorithms do not lend themselves as a robust 
input/output "black-box” [42], they are indeed capable of handling large-scale com- 
putational problems which may be impossible for direct system solvers. It is espe- 
cially desirable to employ the algorithms when one seeks the solutions of large-scale 
numerical systems without resorting to costly computing resources. 

4.2.1 BiCG Algorithm With Preconditioning 

Conjugate gradient (CG) algorithms have been developed for over forty years 
[43,44] and one of the primary applications nowadays is to solve large scale linear 
systems, as aforementioned. It is noteworthy that there exist various versions of 
the CG algorithms, taking advantage of different properties of the matrix such as 
symmetric and sparsity. Also, preconditioning is often used to speed up convergence. 





As sii^gestrcl in the algoritlmi used in tlii> work is a> iollows; 

Cdven 

p, = Ti = b - A X, 

Pi = f, 

For /r = 1.2.3... . 

fl- • Ta- 

Ok - — T 

Pa A ■ pt 

rt+i = rt - qaA • 

ta+1 = ta - «aA*^ • p;,. 

l*A Tit 

Pa+1 = ta + ,^aPa 

Pa+1 = ^k+l + f^kPk 

Xa+1 = Xk + akPk 

where + denotes the complex conjugate and T is the transpose. This version of the 
iterative algorithm is quite general in terms of the system matrix to be solved and 
is usually referred to as the Biconjugate gradient (BiCG) method for unsymmetric 
systems. If the matrix A is symmetric and the initial value is chosen as ri = 
rj, the algorithm can then be shortened to require only one matrix-vector product 
per iteration, since in each step Fa and p;^ are complex conjugate of Fa and p^. 
respectively. 

The ordinary conjugate gradient algorithm can be considered as a special case of 
the BiCG when A is Hermitian (i.e. A = A*^). Again in this case, the algorithm can 
be shortened to have about 50% less computational effort. The CG algorithm is also 
amenable to a straightforward interpretation of its convergence principle. Basically. 
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the algorithm minimizes tlie function 

/(x) = -xA x — b-x 

Hence x oljlained from the CG algorithm after the A- ittualion steps l>e(X)ines the 
solution of the equation 

A/(x) = A - x - b = 0 

The CG type of algorithms did not become competitive until preconditioning 
was introduced to improve the original system condition and significantly reduce the 
convergence rate. One simple preconditioner is the inverse of the original diagonal 
matrix. In our applications this preconditioning has been quite successful and is used 
in conjunction with the BiCG algorithm. 

4.2.2 BiCG— FFT Algorithm For Linear System 

In our work, the linear system of equations is usually large in size, partially full 
and partially sparse. The conjugate gradient (CG) type of algorithms can be used 
to alleviate the memory requirement since the LU decomposition requires excessive 
memory and CPU time. However, the partially dense matrix due to the boundary 
integral equation may still dominate the CPU demand. This is because the method 
of moment (MoM) always leads to a dense system by its nature. Solving the dense 
system in a traditional manner requires 0{N^) order of operations per iteration, 
where N is the boundary integral system dimension. Reduction of the operation 
counts will of course significantly decrease the solution time and this can be accom- 
plished by recasting the BI system onto a few Toeplitz submatrices and making use 
of the fast Fourier transform (FFT) to carry out the matrix-vector products in the 
BiCG solver. 



As (Icsrribccl in tlu‘ next scrlion. the l>uumlaiv intciiial r(|uaiion can i 
a convolutional form if a uniform uricl is api)licd for discn-t i/at k) 1 i. 1 lii" i> mil sui- 

prising since the Green's function is involved in the integration, lu solv the e(iuai lun 
using the CG algorithm, it remains to carry out the convolution at each iKnalion, 
To this end. one may calculate the convolution by taking tin- l ouner Transform of 
two spatial data sequences (arrays) in which the convolution becomes the product of 
the two ’frequence sequences. .An inverse transform of the product \ ields the lesult. 
In contrast to the order of CPU requirement for a matrix-vector product in 

a traditional fashion, the scheme needs 0(.Vlog2A') operation counts if the FFT 
algorithm is employed. The operation reduction is indeed significant and the tech- 
nique has the lowest CPU demand among integral equation solvers, including the 
fast multipole method (FMM) [47], thus is always preferred. 

4.2.3 Convolutional Form of Boundary Integral 


The boundary integral equation is discretized using the structured triangular grid, 
and the relation between the unstructured and structured mesh is described in the 
next section. We recognize that the triangular grid consists of equal right triangles 
as shown in fig. 4.1 and thus involves three different classes of edges (class 1. 2 and 
3). These include the x-directed, y-directed and the diagonal edges, all of which are 
uniformly spaced. For the FFT implementation each class of edges is independently 
numbered in accordance with their geometric location. Specifically, the rth class will 
carry the numbering (m,n) if the edge is the mth along the x direction and the nth 
along the y direction. The indices (m,n) take the v'alues 


m 


0,1,2.... ,M' 


= 0.1. 2..... A” 
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m: 0 1 2 3 M 

Figure 4.1: Structured mesh consists of equal right triangles 


with i = 1 for the y-directed edges, i = 2 for the diagonal edges and 
T-directed edges. Consequent!}', vve find that 


M' 


M - 2 2 = 1 
< M - 1 t = 2 
M - 1 i = 3 



A' - 1 ? = 1 

i\ - 1 t = 2 

N -2 1=3 


•3 for the 


(4.1) 


where M and N denote the numbers of elements along the x and y directions, 
respectively. 

To perform the integrations for the evaluation of the boundary integral matrix 
elements, it is now convenient to rewrite the ba.sis functions (3.14) in terms of the 
new indices (m,n). We readily find that the edge-based basis functions associated 
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(y-nAy)x + (mAx-x)y (.r.y)e5r 1''^ 

0 


otherwise 

where the superscripts refer to the edge class. After the discretization and assemljly 
processes, one obtains a discretized version of the BI system, from winch each entry 


of the boundary matrix-vector product can now be calculated as 

3 

{BI subsystem} = [B] {£*} = ^ 

j=i 

in which (m,n) are the geometric location indices for the zth class observation edges 
whereas ( 777 ',n') are the same for the jth class edges belonging to integration el- 
ements. Thus, the specification of the indices 7, tti and n completely defines the 
^ M of the column resulting after the execution of the boundary 

matrix-vector product. It is readily found that 


AfJ 


sr^ 


Bl 


^ rnn^m' ,n 

m' = 0 n'=0 


• Ei 


(4.o) 


DU 


-2kl SI, ■ SI,,. Go(r, r') dx dy dx' dy' 

IIJL c,(r) Cj(r) l,lj G'o(r. r') dx dy dx' dy 


(4.6) 
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A(/ / = 1 

^/(A.r)- + / = •_> '1.7 I 

A.r / = 3 

More importantly, it can be shown that the BI subsystem ( l.b) exhibits tlie convo- 
lutional property and tlius we can rewrite ( 1..5) as 

3 

[B]{E,) = Y,B'^ * (-J.S) 

j=i 

where the denotes convolution. The proof will be presented in the end of this 

section to ensure the smooth discussion of the remaining procedure. It is now seen 

that the computation of the boundary matrix-vector product can be performed by 
employing the 2-D discrete Fourier transform (DFT), thus avoiding a need to store 
the BI matrix other than those entries which are unique. When the symmetry 
property of invoked, implying 

n'-n) (^*9) 

it is concluded that the total non-redundant entries in the BI matrix are 

3 3 

Ap = - 1) (-I-IO) 

1=1 j=i 

This should be compared to the M' N')^ entries whose storage would normally 

be required if the BI system was not cast in convolutional form and it is notewor- 
thy that Ap IS much smaller in number. To avoid aliasing, it is necessary that 
^{m-m'.n-n') ~ ^ast in a 2-D array w'hich has the usual periodic form, 

and zero padding may also be required to make use of the standard FFT routines. 
Specifically, the malri.x-vector product (4.5) is e.xecuted by using the .MFTxXFT 
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MFT - + 1 < /» < MFT 
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0 otherwise 

with the corresponding field vector given by 


I E^im.n) 0<m<ME 0 < n < -V-' 

(4.1-) 

0 otherwise 

and MFT and NFT must be powers of 2 if a radix 2 FFT algorithm is used. 

In the BiCG-FFT algorithm the BI subsystem vector is symbolically computed 


(BI subsystem} = ^5{DFT-' {DFT(B;J} ■ DFT{£;)}) (4.13) 

1=1 

The presence of the operator 5 indicates the necessary reordering of the 2-D array 
which results after the inverse FFT operation into a single column with the proper 
indexing for addition to the FEM subsystem. It should be remarked that in contrast 
to [40] the integrals (4.6) are evaluated without introducing any approximation. This 
is necessary to preserve the symmetry feature of the global combined system. 

As promised, we now' show (4.8), or the relation 
conclude this section. To simplify the proof, we refer to fig. 4.2 and consider only 
the first integral in (4.6). The same proof can be appied to the second integral m 
(4.6). In addition, with no loss of generality for the proof, the 7 = 1 class edges 
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Figure 4.2: Illustration of two triangles with the corresponding indices to liel]) to 
prove the convolutional property of the boundary integral. 


(v directed) in for the trial function and the ? = 2 class edges (diagonal) in S~ 
for the testing function will be used. To this point, substitutng 5^„(.r.y) in 5+ and 

the first term in (4.6) yields 

= C J J ^ JJ ~ ~ + (x — mAx)(mAx + 2A.r — ,?')] 

Go [(x - x')x + (y - y')y] dxdydx'dy' (4.14) 


where C is a constant coefficient and its detailed form is not of our concern for 
this proof. Note that the integration limits should be set as [mAx.(m + l)Ax] and 
[nAy, (n + l)Ay] for the unprimed coordinates and similarly for the primed ones. 
Therefore, (4.14) will be simplified if the follow'ing transforms are introduced. \ iz. 

X = mAx + ^ x' = mAx -t- 


y = nAy + t] y = nAy + y' 
Indeed, on substituting the transforms, one obtains 


(4.15) 


Intmn.m'n' ~ ^ J J j J 4" Go {.^[(n? — m^)Ax + 

(^-0] + y[(n - n')Ay+ (t) - 7/)]} d^di]d^'dif (4.16) 
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This is the desired relation. From the proof, we can conclude that the convolution 
in (4.8) holds. 

4.3 Mesh Overlay Scheme 

As described above, the BiCG-FFT solver requires uniform aperture griddmg 
so that the BI subsystem can be put in block circulant form. This can be always 
achieved during mesh generation whenever the patches are rectangular in shape oi in 
case of radiators which are placed at some distance (usually small) below the aperture 
as shown in fig. 4.3(a). In the latter case, one might need to add an appropriate 
absorbing material around the edge/corner of the cavity near the aperture to avoid 
possible edge/corner effects, especially when the aperture size is not large enough. 
Fig. 4.3(b) shows the example of this implementation. 

4.3.1 Field Transformations 

However, for circular, triangular, or other non-rectangular patches on the apei- 
ture, it is not natural to construct a uniform mesh using the mesh generator. Typ- 
ically, the aperture mesh is necessary to conform to the patch shape, leading to an 
unstructured free surface grid. In this case, to make use of the efficient, low mem- 
ory BiCG-FFT algorithm, an approach is proposed to overlay on the unstructured 
aperture grid another coincident structured grid as shown m fig. 4.4. The boundary 
integral subsystem is then constructed by using the overlaid uniform grid whose edge 



uniform grid 



(a) 

Circular Patch Recessed Underneath the Ground Plane 
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(b) 

Figure 4.3: Printed circular patch antenna is modeled using the recessed scheme to 
incorporate the BiCG-FFT algorithm, (a) Illustration of the configura- 
tion; (b) Comparisons of the BiCG-FFT result with that of the ordinary 
FE-BI technique presented in chapter 3. 








Figure 4.4: Overlay of a structured triangular aperture mesh over an unstructured 
mesh, shown here to conform to a circular patch. 

fields can be related to those on the unstructured grid via two sparse transformation 
matrices. That is, it is necessary to append to the system (3.11) the relations 

{Es}^ = [^f] 

= [TB]{Esh 

w'here the subscripts u and nu refer to the field coefficients of the uniform and 
non-uniform aperture grids, respectively. Also, [7>] and [Tg] refer to the forward 
and backward transformation matrices, respectively, with A'u and A'„u denoting the 
numbers of the uniform and non-uniform mesh edges on the cavity aperture. 

To derive the elements of [Tg], we begin with the expansion (3.14) and enforce it 
at three points on each edge belonging to the uniform grid. We convenient!}- place 
these three points at the center and ends of the edge (see fig. 4.5). 

Given the fields at these points, we can interpolate the field along the {m.n ) edge 






overlaid unstructured 
mesh triangles 
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Figure 4.5: Illustration of the parameters and geometry used in constructing the 
transformation matrix elements between the structured and unstructured 
mesh. 


of the uniform grid using the weighted average 
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(4.19) 


in which denotes the unit vector along x,y or the diagonal, depending on the class 
of edge being considered. The quantities represent the fields in the non-uniform 
grid triangles with the superscript k being a sum variable in case Tend] . fendj or Tmid 
specify a point shared by more than one triangle. Obviously. Aend, < A'mid and .Vend^ 
denote the number of non-uniform grid triangles sharing the node at Fendj • Tmid and 
fendj, respectively, and will typically be equal to unity. 

.After assembling (4.19) into (4.18) we find that the elements of the forward 
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in which 


1 J = j> 

0 otherwise 


and the global indices i and j correspond to the ?th uniform grid edge and the jt h non- 
uniform grid edge. The subscript j( is the global index used in numbering the non- 
uniform grid edges, whereas the subscript f (= 1, 2 or 3) is the local edge index used 
in the definition of the basis functions Sf. We remark that the explicit computation 
and storage of the transformation matrix elements results in a substantial increase in 
efficiency because it avoids the usual assembly process during each iteration step and 
that the proposed overlay scheme allows the analysis of large non-rectangular patch 
arrays because storage of fully populated BI system matrix is avoided. The user 
needs to only provide an additional data file which flags the uniform grid edges lying 
on a PEC element and this is an important user-oriented feature of the formulation. 
Following the same procedure, we can obtain the expression for the entries of the 
backward transformation matrix, ft should be noted that assuming each uniform 
grid edge traverses three or less non-uniform grid triangles, the non-zero entries in 
each row' of [Tf] will be 9 or less. However, they can reach a maximum of 18 if the 
midpoint and endpoints reside on an edge of the non-uniform grid. The maximum 
non-zero entries in each row of [Tb] will be 15. but the typical number is much less. 


4.4 Results 


Figure 4.6 shows a cavily-backcd 2 x 2 patch array, where ea< h patc h is a riiihi 
angled triangle. Since this geometry is ada[)table to a uniform mesh with riiilii 
angled triangles, it is used to verify our proposed FE-BI techniciue incorporated 
with the BiCG-FFT system solver. The developed FE BI code with the BiCCJ I'FT 
is first compared with the original version of the hybrid FE-BI techniciue dc'scril)c‘d 
in chapter -3. As shown in fig. 4.7. the monostatic radar cross section (RCS) ])attern 
over the space 0 < 6* < 90° at the o = 0 plane agrees very well with that compiitc'd 
using the regular BiCG solver. 

It is also informative to compare the scattered fields by the same cavity-backed 
structure without the patch array to find the contribution of the array to scattering. 
Figure 4.8 shows the monostatic RCS patterns by the aperture with the absence of 
the patch array. Again, the computations were obtained using both the regular BiCG 
and BiCG-FFT v'ersions of the FE-Bl methods. As can be seen, the level of the 
scattered field at the normal incidence reaches above zero in dB/A^ with the presence 
of the patch array, whereas the scattering by the aperture at the same incidence is 
about 23 dB/A^ lower. 

To varify the overlaying scheme for nonrectangular geometry, we evaluated a 
bistatic RCS scattering as shown in fig. 4.9 by a cavity-backed circular patch an- 
tenna. In this case, the dielectric fillings of = 4 and = 10 inside the cavity were 
used, respectively, and the results obtained using the regular BiCG FE-BI method 
are compared with those computed using the BiCG-FFT with the overlaying trans- 
formations. It is observed that the agreement is quite satisfactory in scattering 
analysis. 



Patch Array 



Figure 4.6: Illustration of the cavity-backed triangular patch array. 

For radiation analysis (e.g. input impedance) where an accurate prediction of 
the near-zone fields may be required, the accuracy of the overlaying sclieme can be 
significantly enhanced by considering the trial-testing element's interactions in the 
boundary integral. Specifically, it is suggested to separate the interactions between 
the closed-region elements from the far zone weak couplings. The strong close- 
region couplings are treated using the normal method of moments, whereas the weak 
coulings are computed using the fast algorithm. This approach has been reported 
(see e.g. [47]), and once combined with the overlaying scheme, it can be used to 
control the accuracy of the FE-BI technique in an adaptive manner. 





A 

Figure 4.7: Comparisons of the monostatic radar cross section scattering by a 2 x 2 
triangular patch array shown in fig. 4.6. The reults were computed using 
the regular BiCG FE-BI technique described in chapter .3 and using the 
BiCG-FFT proposed in this chapter. 
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Figure 4.8: Comparisons of the monostatic radar cross section scattering by an empty 
aperture with the same cavity size and dielectric filling (er=l) as the 
structure shown in fig. 4.6. Again, the results were calculated using 
the regular BiCG FE-BI technique described in chapter 3 and using the 
BiCG-FFT proposed in this chapter. 
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Figure 4.9; Bistatic RCS scattering by a crcular patch antenna modeled using the 
regular BiCG FE-BI and the BiCG-FFT algorithm with overlaying 
transform. 




CHAPTER V 


Efficient Finite Element Subsystem — II 
5.1 Introduction 

As demonstrated in the previous two chapters (also reported in [13.29.48]). a 
hybrid finite element-boundary integral technique [48.13] can be employed for char- 
acterizing conformal antennas of arbitrary shape [29]. Indeed, planar/non-planar. 
rectangular/non-rectangular designs, ring slot or spiral slot antennas with probe, 
coax cable or microstrip line feeds can be simulated with this formulation. This 
is because of the geometrical adaptability of tetrahedral elements used for the im- 
plementation. However, in practice, certain configurations require extremely high 
sampling rates due to the presence of fine geometrical details. Among them are a 
variety of slot antennas (spirals, rings, slot spirals, cross slots, log-periodic slots, 
etc.), where the slot width is much smaller than the other dimensions (cavity diam- 
eter or inter-distance of slots). In these cases, the mesh is extremely dense (with 
over 50. 100 or even higher samples per wavelength), whereas typical discretizations 
involve only 10-20 elements per wavelength. This dense sampling rate is especially 
severe for 3-D tetrahedral meshes, w'here the geometrical details usually distort the 
tetrahedrals. The numerical system assembled from this type of mesh often leads to 
a large system condition due to the degraded mesh quality. .Also, mesh generation is 
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Figure 5.1: Geometry of cavity-backed microstrip antennas 

tedious and the solution CPU time is unacceptably large. 

In this chapter, we propose a finite element-boundary integral formulation using 
edge-based triangular prism elements. It can be shown that this element choice is 
ideally suited for planar antenna configurations containing spirals, circular and trian- 
gular slots. .Among many advantages of the prismatic elements, the most important 
is the simplicity of mesh generation. Also, much smaller number of unknowns is 
required for an accurate and efficient modeling of complex geometries. Below, we 
begin by first outlining the finite element-boundary integral (FE-BI) formulation for 
slot antenna modeling. A new, physically meaningful, set of edge-based functions 
for prisms is then presented to generate the discrete system of equations. Finally, 
the applications of the proposed hybrid FE-BI method to various antenna radiation 
and scattering problems are given to conclude the chapter. 

5.2 Hybrid FE-BI Formulation 

Consider the cavity-backed slot antenna shown in fig. 5.1 where the cavity is 
recessed in a ground plane. To solve for the E-field inside and on the aperture of the 


cavity, a standarci aiJjJioacli i> lc» (‘.\t remizc t lu* fiinct iunal i M > wtiich. iui i .uli.ii iuii 
ami sraneriiig jjrol)lcnis. may be generalized a> 

F(E) = \ xE)-ji:' -(V ^ E)- k-Ef. E}J\ 


+ 


E • (jA-uZoJ' + V X j?:' . M‘) 


<l\ 


+ jkoZo [ I E ■ ( H X n ) (IS 
J Jso+Sf 


(ml) 


where and JI^ denote the relative tensor constitutive parameters of the ca\itv 
medium. Zo and kg are the free space impedance and j)topagation constant, respec- 
tively, 5’o represents the aperture (or slots) excluding the metallic portions and Sj 
denotes the junction opening to the guided feeding structures. .Also, f j is the volume 
occupied by the source(s) (J' or M') and H is the corresponding magnetic field on 
5'o and Sj whose outer normal is given by n. As before, the explicit knowledge of H 
in (5.1) is required over the surface Sq and Sj (also referred to as mesh truncation 
surfaces) for a unique solution of E. Specifically, the magnetic field H over .S’o may be 
replaced in terms of E via a boundary integral (BI) or absorbing boundary condition 
(ABC), whereas H on 5/ is determined on the basis of the given feeding structure. 
This version of the functional as compared to (3.1) allows an easy inclusion of various 
feed models, such as aperture coupled slot, coax cable, etc. (see chapter 6 for details). 
In this chapter since we concentrate on improving the FEM efficiency, therefore the 
boundary integral method will be employed as in chapter 3. It will be seen that this 
implementation indeed meets the accuracy need without extra CPU burden. In the 
context of the FE-BI, H is represented by the integral 

H = ojkoYo J I Go(r,r') • (f x E(r')) ds’ 


: 5 . 2 ] 


where G is the electric dyadic Green's function of the first kind such that f? x G = 0 ie 


satisfied on the (planar, spherical or cylindrical) metallic platform (refer to chapter 1 ), 
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with r and r' being the observation and integration points. resp(>r1 ivel\ . ami I = 
H- + j/y + ff is the unit dyad. In connection with our problem, i.e. that of a cavity 
recessed in a ground plane. is equal to the sum of the incident and reflect ('d 

fields for scattering computations, or zero for antenna parameter evaluations. 

To discretize the functional (5.1). we choose to subdivide the volume region using 
prismatic elements as shown in fig. 5.2 and fig. 5.3. The field in each of the prisms 
can be approximated using the linear edge-based expansion [49-51] 




( 5 . 4 : 


J=l 


where [V]e = [{V;}, {V;}, {V;}], and {E^} = {£,% E ^. . . , E^V . The vectors u : 

x,y,z, are of dimension m = 9 and they simply represent the x,y,z components of 
Vj associated with the jth edge of the eth element. Since are chosen to be edge- 
based functions, the unknown coefficients represent the average field along the 
jth edge of the eth element. A corresponding representation for the aperture fields 


E(r) = J]£fS*(r) = [5]r{£*}, (5-5) 

z = l 

where [5]. = [S'x,5j. and V(r) reduces to S"(r) when the position vector is on the 
slot. 

To generate the discrete system for Ej, (5.4) and (5.5) are substituted into (5.1) 
and subsequently f(E) is differentiated with respect to each unknown Ej. With the 
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where the sums are over the total number of volume or surface elenumts. In this, tin 
matrix elements are given by 

V X V, - A-;^V E ■ V| (Ir 


A7 = J^oZo J V, ■ 


3. i 


j- + V X ■ V X M' f/f 


(0.8) 


= - JJ JJ 2A^S*(r) • S;(r')G'o(r,r')cfs</s' 

+2 [V X S^(r)],[V' X S^^{r')lGo(r.r')dsds' (5.9) 

17 = 2jkoZojj^S:-{H'x=)ds (5.10) 

where the subscript c in denotes to take the z component. It is noted that L - 
is removed in case of radiation problems and that the same holds for A, when the 
scattering problem is considered. 


5.3 Edge-Based Prismatic Elements 


Consider the right angled prism shown in fig. 5.3 w'hose vertical (z-directed) sides 
are parallel (right-angled prism). We now design two geometric quantities as 

hi = e, ■ z X {r - r,), S- = -^h, (5.11) 

where r, denotes the location of the fth node, e. is the unit vector along the fth 
triangular edge opposite to the fth node, denotes the length of this edge and r 
is any position vector terminated inside the triangle. One way to obtain an edge- 
based field representation for the prism is to utilize the nodal basis functions [52] 
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and t]i(*n a[j[j|\- tlic prorcdun' disriissf'd in [ id. od. '» l'. llouvvrr. an allcmal i\c .md 
more j)liysically meaningful approach can l»e ein|»loyed fur llie < <.)n>i nicl ion of ilic 
edge elements. Referring to fig. -j.2. it is evidtmt that if r is in the x > plane, then 
5'ff in (5.11) gives the area of another triangle 12'3' sucli tliat the lengths of edges 
joining the nodes 2 — 3 and 2' — i' are equal. With this definition of r. wc define a 
vector a vector 


S- 

where i, ■ S, = That is. the vector component along c, has a magnitude which 
is equal to the ratio of the areas of the triangle 12'3' to that of 123. We observe 
that (5.12) is simply the edge-based expansion for the triangular elements [32] and 
is the appropriate expansion to be used in (5.5). The corresponding volumetric basis 
functions can be obtained by inspection, viz. 

V, . = 1.2.3 

j = 4,5,6 (5.13) 

k = k = 7,8,9 

where Ct is the triangle simplex coordinate associated with the A’th prism vertex at 
(^k-,yk)- As illustrated in fig. 5.3, Zc and h = A- represent the offset coordinate and 
the prism height, respectively. When (5.13) are substituted into (5.7), the resulting 
integrals can be evaluated in closed form as given in the Appendix C. However, 
the integrals resulting from the substitution of (5.12) into (5.9) must be carried out 
numerically, except the self-cells which can be performed analytically as discussed 
by Wilton [55]. 



5.4 Applications 


Thin slot antenna structures ha\e been treated usinu the al.ovc' formulated 1 1. HI 
technique and certain modeling results will he presented m this section to (Uunon- 
strate the validity and capability of the approacli. 

Raciiation and scattering by an Annular Slot; 

To evaluate the accuracy and efficiency of the prismatic mesh and the afore- 
mentioned implementation, we first consider the analysis of the narrow annular slot 
(0.75cm wide) shown in fig. 5.4. The slot is backed by a metallic circular cavity 24.7 
cm in diameter and 3 cm deep. The FE-BI method is quite attractive for this ge- 
ometry because the slot is very narrow and most of the computational requirements 
are shifted on the finite element portion of the system. The calculation shown in 
fig. 5.5 were carried out using the prismatic and tetrahedral elements [29], .^s seen, 
they overlay each other. However, only 1024 prisms were needed for modeling the 
cavity, whereas the number of the tetrahedral elements for this homogeneously filled 
cavity were 2898 for acceptable element distortion. If a multi-layered structure was 
considered, or a similar system condition was used as a criterion for mesh generation, 
then much more tetrahedrals than prisms would be needed for modeling such a struc- 
ture. Moreover, the prismatic mesh is trivially generated given the slot outline. In 
contrast, substantial time investment is required for generating and post-processing 


the tetrahedral mesh. 

Frequency Selective Surfaces: 

Frequency selective surfaces (FSS) structures [56,57] are arrays of tightly packed 
periodic elements which are typically sandwiched between dielectric layers. The 
periodic elements may be of printed form or slot configurations designed to lesonate 




Figure 5.4: Geometry of the annular slot antenna backed by a cavity 24.7 cm in 
diameter and 3 cm deep 

at specific frequencies. As such, they are penetrable around the element resonances 
and become completely reflecting at other frequencies. To meet bandwidth design 
specifications, stacked element arrays may be used in conjunction with dielectric 
layer loading. 

Here we shall consider the analysis of FSS structures (with slot elements) via 
the FE-BI method. Because of the fine geometrical detail associated with the FSS 
surface, the finite element method has yet to be applied for the characterization of 
FSS structures, but use of prismatic elements makes this a much easier task. Of 
particular interest in FSS design is the determination of the transmission coefficient 
as a function of frequency, and since the array is periodic, it suffices to consider a 
single cell of the FSS. For computing the transmission coefficient T, the periodic cell 
is placed in a cavity as shown in fig. 5.6 and the structure is excited by a plane wave 
impinging at normal incidence. Assuming that near resonance the wave transmitted 
through the FSS screen will retain its TEM character, the transmission line concept 
can be used to find the scattered field 

aT^ 

1 — qR 


E* = 







Figure 5.5: Scattering: Bistatic (co-pol) RCS patterns computed using the tetrahe- 
dral FE-BI code and the prismatic FE-BI code. The normally incident 
plane wave is polarized along the </> = 0 plane and the observation cut 
is perpendicular to that plane. Radiation: X-pol and Co-pol radiation 
patterns in the 0 = 0 plane from the annular slot antenna shown in 
fig. 5.4. The solid lines are computed using the tetrahedral FE-BI code 
whereas the dotted lines are computed using the prismatic FE-BI code. 
The excitation probe is placed at the point (y=0) marked in fig. 5.4. 
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Figure 5.6: Illustration of the setup for computing the FSS transmission coefficient 
Upper figure: periodic element (top view); Lower figure: periodic element 
in cavity (cross-sectional view) 
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where* 7 is the transmission roe'fheient of the* I SS. f\ = \ — I ami i> i hr rt'llrr- 
lion coefficient associated with the ca\'it>‘ base. To reduce the multiph* inttu cu i luii'- 
within the cavity, it is appropriate to terminate the cavity with soiiu' absorber, ihii'' 
reducing lire value of o to less than 0.1. Since H is also small near resonamc. a good 
approximation for T is 

C = ioi°e 

and upon considering the next higher order cavity interactions, we have 

r,B~rj°> + ioiog[i-a(i-r<°>)]. 

A more direct and traditional computation of Tjb would involve the placement of 
the FSS element in a thick slot [58]. However, this requires enforcement of the 
boundary integral over the entire lower surface of the slot, leading to a much more 
computationally intensive implementation. 

The above FSS modeling approach was applied for a characterization of single 
layer and multi-layer FSS structures. In both cases, the periodic element was a slot 
configuration. The geometry of the single layer periodic element is shown in fig. 5.6 
and consists of a planar slot array on a dielectric layer 0.0762 cm thick and having 
Cr = 4.5. The FE-BI calculation using prismatic elements is given in fig. 5.7. Clearly, 
our calculations are in good agreement with the measurements and data based on 
the more traditional method of moments [59,60]. 

The geometry of the multilayer radome considered in our study is given in fig. 5.8. 
The total thickness of the FSS was 6.3072 cm and is comprised of two slot arrays (of 
the same geometry) sandwiched within the dielectric layers. For modeling purpose, 
a 1.54cm thick absorber is placed below the FSS as shown in fig 5.8. From the 
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Figure 5.7: Calculations and comparisons of transmission through the FSS structure 
shown in fig. 5.6 

calculated results, it is seen that the results generated by the FE-BI method are in 
good agreement with the measurements. 

Radiation Property study of Conformal Slot Spiral Antenna: 

Consider a typical Archemidean slot-spiral antenna shown in fig. 5.9. This an- 
tenna is built on a double-sided PCB with its two arms following the expression: 
r = Q0 /3, where a = 0.1333cm and (3 = 2.8595cm. One arm can be determined 
from the other by rotating 180° counterclockwisely. It is noted that this structure 
differs from the conventional design in that the central portion of the spiral is not 
fabricated. The reasoning for it relies on the facts that the antenna is designed with 
a bandwidth less than 30%, and that the central portion usually requires a careful 
fabrication because of the geometric details, and still that the central space may be 
used for possibly complex feed network. One of our goal is to study the effect of this 
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spiral shape on its [x'rforniancc. 

A henclirnark test tnodel is designated to operate from IIS Mil/ to 1 *>7 Mil/, 
to replace the con\ ent ional protruding hlade antenna. I'lie size h(.)\ve\cr is niii(h 
compact with its conformality property. Our simulation model is scah-d h\ I s tu 
operate at 944MHz to 12o6MHz with the center frequency llOOMHz. The values of 
o and 3 above were determined based on this frequency band and also the iiuiiiIk'i 
of turns (4.5). The cavity is filled with a dielectric slab (<r = 2.2) of 0.3 cm depth, 
corresponding to approximately 0.011 free space wavelength at the center frequency. 
The antenna's directivity is analyzed from the radiated pattern at lower, center and 
higher frequencies and the results are tabulated in Table 5.1. 

Figure 5.10 - 5.12 show the radiation patterns for frequency 944. 1100 and 1256 
MHz, respectively. The Eg and at the principle plane (p = 90° are plotted. It 
is understood that when the frequency varies, the active region travels along the 
slot spiral. Thus the principle plane may not be coincident with the E-plane. (In 
fact, the E-plane is not clearly defined in this case.) The optimum axial ratio for 
the three cases are tabulated also in Table 5.1, and it shows that the spiral shape 
design really plays an important role to insure a good quality radiation pattern. At 
both center and lower frequencies, less than 3 dB AR has been achieved. When the 
frequency increases, the active region moves inwards to the center and becomes closer 
to the feeds where the EM fields exhibit a comparatively strong profile. The radiated 
pattern therefore is most likely affected and this explains why the AR increases at 
the high frequency. The AR deterioration can be avoided by adding a couple of spiral 
turns inside. It is seen, nevertheless, that a CP mode can be achieved within the 
entire designated bandwidth and with a wide azimuthal angle (as wide as 60° in the 
optimum case). In practice, we notice that absorbing materials may be needed to 
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regulate the magnet ir currents at the Ix-oinninc or eml" of the 'lot 
when the iiuinl>er of turns is mininiiz<‘cl. 
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iTequeiicy ((IHz) 

j U.TO 

1.100^ 

1.2')(. ; 

Ciain (dB) 

j 7.-J2 

(i.bb 

5.22 ' 

.-\xial Ratio (dB) 

2.7 

1.0 
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Table 5.1: Comparisons of gain and axial ratio at different oi)erating fr(Hiu('iu ies 


5.5 Concluding Remarks 

A hybrid finite element-boundary integral (FE-BI) formulation was presented 
for modeling narrow slots in metal backed cavities. Prismatic elements were used 
in connection with the FE-BI implementation, and in contrast to the tetrahedral 
elements, these offer several advantages. Among them, low sampling rates are needed 
for generating meshes and the mesh generation process is substantially simplified. 
Other advantages of the prismatic elements over the tetrahedral elements include 
better system conditions and faster pre/post data processing. 

The explicit expressions for FE-BI implementation of prismatic elements were 
tabulated and numerical results for slot antennas and frequency selective surfaces 
were presented to demonstrate the validity and capability of the technique. 
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Radiation Pattern at f=0.944GHz (lower end of frequency range). It 
can be seen that the axial ratio of the pattern becomes larger comjiared 
to that at the center frequency, but still remains within 3dB for a wide 
angle range. This indicates that the number of the outer turns in the 
spiral contour design is most likely sufficient. 


frequency* 1 .256 GHz 

- - 0 - --- 



solid line: E_phi. dashed line: E_theta 


Radiation Pattern at f=1.256GHz (higher end of frequency range). It 
can be seen that the axial ratio of the pattern is deteriorated compared 
to those at the center frequency and lower frequency. This certainly 
shows that the number of inner loops still needs to be increased to 
insure a good quality pattern. 




CHAPTER VI 


Antenna Feed Modeling 


For scattering problems where the plane wave incidence is usually considered as 
the ‘source", the right-hand-side excitation has been explicitly expressed in (3.7) and 
(5.10) and will not be discussed further. However, for antenna impedance evaluations, 
we have proposed and reported several feeding schemes [61] associated with various 
practical feed designs for microstrip antennas. Some of these are discussed below. 

6.1 Probe Feed 

6.1.1 Simple Probe Feed 

For thin substrates the coaxial cable feed may be simplified as a thin current 
filament of length / carrying an electric current II. Since this filament is located 
inside the cavity, the first term of the integral in (3.7) needs to be considered for this 
model. Specifically, the zth (global) entry of the excitation vector A'j becomes 

A', = jkoZoI i ■ V,(r), i = ji, j 2 , Jm 

where r is the location of the filament, m is the number of (non-metallic) element 
edges and jm is the global edge numbering index. In general, m such entries are 
associated with m element edges, and thus the index i goes from ji up to . This 


98 



cxprossioii can 
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coincident with tlie current filament. 

6.1.2 Voltage Gap Feed 
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This excitation is also referred to as a gap gt in valor nud amounts to spenfvinu, a 
priori the electric voltage \ ' across the opening of the coax cable or any othei gap. 
Since V = E ■ d. where d is a vector whose magnitude is the gap width, and E the 

electric field across the gap. we have that E, = " ^ 

?lh edge is parallel to d. Numerically, this gap voltage model can be realized by first 
setting the diagonal term .4,, equal to unity and the off-diagonal teims .4,j [i ^ j) 
to zero. For the right-hand-side vector, only the entry corresponding to the nh 
(global) edge across the gap is specified and set equal to the value E, whereas all 
other entries associated with edges not in the gap are set to zero. 


6.2 Aperture-coupled Microstrip Model 


As shown in fig. 6.1, when the microstrip antenna is fed with a microstrip line 
network underneath the ground plane (cavity’s base) via a coupling aperture, special 
treatment of the feed structure must be considered in the FEM formulation. This 
is because the microstrip line is usually designed to have different size and shape as 
compared to the cavity’s geometries. Hence, the conventional simulation of treating 
the entire 3-D domain using a single type of elements is not efficient or appropriate 

for this feed. 

Referring to fig. 6.1, it is appropriate to separate the computational domains 
because of the small element size required in modeling the guided feed structure. 
One difficulty encountered when this decomposition is implemented is how to model 
the coupling through the aperture. As an example let us consider a rectangular 
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Figure 6.1: Cross-section of an aperture coupled patch antenna, showing the cavity 
region I and the microstrip line region II for two different FEM compu- 
tation domains. 


(a) 


(b) (c) 

Figure 6.2: slot and its discretization (a) slot aperture; (b) typical mesh from cavity 
region; (c) uniform mesh from microstrip line region. 

edge-based field e.xpansions, the meshes across the common area (coupling aperture) 
are different, and this causes difficulty in enforcing field continuity across the slot 
aperture. However, since the aperture is very narrow, a ‘static’ field distribution 
may be assumed at any given frequency. Therefore, the potential concept may be 
again applied to relate the fields below and above the aperture. Specifically, the 




-potriil ial rout iiniity roiiditio!i is (“iitorr<‘tl. ami to |>ror<vtl u> du so. let u> limt 
classify the slot edges as follows 

Tetrahedral Mesh (Cavity Region I): 

• j = 1. 2. 3. ... vertical edges 

• £”^2 j _ 1.2.3.... diagonal edges 

Brick Mesh (Feed Region II): 

• 7 = 1.2.3,... vertical edges 

Then the ‘equi-potential' continuity condition requires that 

= CjE; (6.1) 

£f = + <,+.£'+.) (6-2) 

in which 

whereas t and d are the lengths of the vertical and diagonal edges, respectively. 
That is, t is simply the width of the narrow rectangular aperture. The coefficient Cj 
is equal to ±1 depending on the sign conventions associated with the meshes above 
and below the coupling aperture. 

The connectivity scheme for entirely different computational domains may be 
extended by generalizing this concept. It is apparent that this approach makes 
the FEM implementation straightforward for different geometry/size domains that 
would be significantly inefficient if only one type of elements were used for modeling 
the structure. In addition, the technique ensures a good system condition since the 
number of distorted elements in the mesh are minimized. 



6.3 Coax Cable Feed 


6.3.1 Motivation 

The coax cable is widely used as a fe('diiig strucUire for microst ii|> oi ia\it\- 
backed patch antennas because of its simplicity and low spurious radiation. Indeed, 
abundant literature exists on the theoretical and ex|)eriniental investigation of coax 
cable feeds [62-64]. Most of these papers present integral e{|uation techniques in 
conjunction with the pertinent Green's functions. However, the Green's function is 
only available for a certain class of geometries, and this limits the application of 
the integral techniques to those geometry designs. .Mso. the formulation must be 
modified and recoded for different antenna configurations corresponding to Green's 
function variations. To avoid the complexity of the Green's function, we recentlv 
proposed a hybrid finite element - boundary integral approach [29] which is described 
in chapter 3 and 4. For antenna radiation, it is observed that a simple probe model 
with a constant current along the inner conductor linking the grounded base to the 
antenna element is straightforward and efficient. But the probe feed is only valid 
for thin substrates and this is consistent with the Moment Method (MM) results. 
To model an electrically thick substrate, in this section a more sophisticated feed 
modeling scheme is proposed in the context of the finite element method (FEM) 
using linear edge-based tetrahedral elements. The formulation of the entire hybrid 
numerical system will be first described in the presence of the necessary functional 
term for feedline. The proposed feed model is then presented on the basis of a TEM 
mode excitation. Model improvements are also discussed for the case when the TEM 
assumption at the cavity-cable junction does not hold. 



6.3.2 Hybrid FE-BI System 
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The functional peilinent to tlie radiation by a cavit\ l)ack<-d antenna with a t uax 
cable feed (as shown in fig. 6.:f) is given by 

F(E) = i |(V X E) ■ ^(V X E) - A'i^t.E ■ e| dc 

-2kl Jj^(E X £) •{/£(£ X i) • ^'u(r.r')d.s'| ,l> 


-jf^oZojJ {ExH)-=dS. 


((i.:r 


where V refers to the cavity volume and the surface S encompasses the cavity apertuta' 
excluding the portion occupied by the metallic antenna elements: ty and denote, 
respectively, the relative permittivity and permeability; Aq is the free space wave 
number, Zq is the free space intrinsic wave impedance, I is the unit dyad, and 
G’o(r,r') is the free space Green’s function wdth r and r' denoting the observation 
and integration points; the surface C is the cross section of the coax cable at the 
cavity-cable junction. 

Following the standard discretization procedure [29], we obtain a system of equa- 
tions of the form 

A> A'se 


E + E {isa(E;}} + e = <>. 


(6.4) 


e=l e€^ 

w'here the explicit expressions for A,j and may be found in [29] and the functional 
term Fc is the surface integral on C in (6.3). 

6.3.3 Proposed Coax Feed Model 


To proceed with the evaluation of 

Fc = -jkoZoll^{Ex}i)-US. 


(6.5) 


a houiKlar\ constraint relating E to H is needed. 'It) this end. uc .i--iime <i 1 l.M 
mode on (' and conse(nient ly the fields witliin the (a\ity max l>e e.xpres'cd a- I'.'e 
hg. (i.d ) 


E = ^ ^(l + r)-r. H = — (l-l)-o. 


2~y/t 


■)- 




where trc the relative permittivity inside the coax cable; F deiu)l('s the n'flectit)n 
coefhcient measured at r = 0 and /o is the given input current source at the same 
location. Also, (r.o.:) are the polar coordinates of a point in the cable with the 
center at r = 0. To simplify the analysis, we introduce the quantities 
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Hence, 


and from (6.7) it follows 




( 6 . 8 ) 


Ao=-^e„+^, (6.9) 

Zo 7T 

which is the desired constraint at the cable junction in terms of the new quantities 
ho and cq. Note that eo and are field coefficients as new unknowns in place of the 
fields E and H, and it is therefore appropriate to rewrite Fc in terms of these new 
coefficients. To do so. we substitute (6.7) and (6.9) into (6.5) and upon making use 
of the axisymmetric field property we obtain 


Fc = -2'!TjkoZoeoh"''ln{-), (6.10) 

a 

where a and b are the radii of the inner and outer cable conductors. The superscript 
src stands to indicate that ho is treated as a source term in the extremization of the 


functional. 
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\\c choose the liiK'ar (‘dge-basecl tetraliedral el<‘uients to disereii/t- the 
and tin* corres|)ondiiig mesh on t!i(' cross section (' is shown in iig. In thi'- 

formulation, the field across the //' edge. p=1.2 \( ( .\ r is the immlxT of ( a\ ily 

mesh edges on C'j- is set to a constant as dictated by the linear etlge-based expansioti 
function inside the cavity. However, the cable TEM modal fields ((>.(>) behave as 
]/r and this modeling inconsistency makes it difficult to apply the tangetitial field 
continuity condition at the cable junction ( i.e. over the aperture (’). To overcome 
these difficulties, we can relate the fields across the cable junction by recognizing that 
the potential difference between the inner and outer conductors must be the same 
as computed by the fields of the cavity or those in the cable region. Specifically, 
if the edge of the cavity region borders and is also across the coax cable, from 
(6.6)-(6.8) it then follows that the appropriate equi-potential condition is 

AV' = Ei{b — a) - eo/n-, i = Np{p — 1,2, ....Ac)- (6.11) 

a 


where AV ' denotes the potential difference between the inner and outer surface of the 
cable. This condition simply provides a relation between the constant cavity edge 
field and the coax cable modal field. When used into the functional Fc, it introduces 
the excitation into the hybrid finite element system without a need to extend the 
mesh inside the cable or to employ a fictitious current probe. It remains to rewrite 
Fc in terms of F,, i.e. the field value of the edges bordering the cable and to do so 
we substitute (6.9) and (6.11) into (6.10). Then upon taking into account all A'c 


cavity mesh edges on the cable junction, we obtain 


^ ^ ^ .1 N j ^0 y/tFcb-a 

Fc = -^3k,Zo{b 


F. £f.. 


( 6 . 12 ) 


In this expression, rather than representing the functional Fc in terms of a single edge 


field, we made use of the average field across the cable as computed by the totality 
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of tlif* equal element fields on the cable's aperture i because u 
pro]>erty. all element^ fields at the cable's ajierture ar<- e<]ual i. 
the curly lirackets of (fi.ld). with the sujierscript >rc. functions 
extremization process. Hence, the extreniization of (O.ld! yields 
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/<77 b — a 


Zo In- 
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L\E, - 7 = -Vp(/>= 1.2 \'c ). 


E, 


f tli<' axisyiiinu't lie 
I he lactoi ln^lde 
as a source in i he 


(b.ld) 


where 


.1 ^{b-af 

a 

(6.14 

- j'^koZo(b — g)/o. 

(6.15 


We observe that the ‘constant cavity field’ along each mesh edge at the cable junc- 
tion is just a fictitious field representation and its meaningful physical interpretation 
is governed by the equi-potential constraint (6.11). To proceed, we assemble the 
FEM system together with (6.13). Specifically, each V, is added to the Nc diagonal 
entries of the finite element matrix which is associated with the Nc edges bordering 
the coax cable. Also, the excitation column of the hybrid system is nullified ev^erv- 
where except for the Nc entries which are set to Vj. Once the hybrid FE-BI system 
is solved [29], the input admittance at z=0 is calculated from 


y 


tn 


TT <f H • r rd4> 


2/o 1 

eoln^ Zc 


(6.16) 


where Zc is the characteristic impedance of the coax cable. 

In the above feed model we assumed the presence of only the dominant (TEM) 
mode at the cavity-cable junction, an assumption which may not be suitable for 
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certain a])plicat ions, lo (>\('iconie this limitation, one approai h i^ to e.xten.i the 
nu-sh jsav. a distance d) intcj the calde. The e(|ui-iH)tential condition udi then l>e 
applied at z=-d. where all higher order mod<‘s vanish. 1 his M'heme recpiin'^ a nioK- 
suitable expansion for the fields in the -d < z <Q section to avoid the comi)li< at ion 
of extending the tetrahedral mesh into the cable and. thus, retain the elliciem vof tin- 
equi-potential feed model. Since in most cases the antenna is operated in a frecpn'iicy 
range far below the cut-off of the first higher order mode of the coax cable, tlu- fi('ld 
distribution near the junction C will still be dominated by tlie fundamental I KM 
mode. With this understanding, a possible suitable expansion for the field in the 
coax cable (using shell elements rather than tetrahedrals) is 


9 • 


(6.17 


where q=r, 6 or z. i= 1.2,3 or 4 and N^(r,(p,r) is the shape function for each of the 
12 edges (3 directionsx 4 edges per direction). They are given by 




(6.18) 


with qai<ib S'lid 9c representing r,(f) and 2 in cyclic rotation and correspondingly 
Qa-Qb and qc represent the parameters r, ^ and c. Also, i denotes the edge number 
along each coordinate, and Aqa is the width of the edge along the direction. The 
correspondence between the edge numbers and the node pairs for each coordinate! r. o 
or z) is given in Table 6.1 along with the definition of the tilded parameters m (6.1b). 
When an axisymmetric field property is assumed, the numerator of the expansion in 
(6.17) reduces to the standard brick element format for the radial and z components, 
independent of the 4> variable. Note also that the particular property of this expan- 
sion is the introduction of the 1/r factor, simulating the coaxial cable mode. The 
accuracy of (6.17) is demonstrated in fig. 6.5. where we show that only 2-3 elements 
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Table G.l: The correspondence between the edge numbers and the node ])airs for each 
coordinate(r. o or z) along with the definition of the tilded parameters in 
(6.18). 


are needed along the radial direction for the accurate prediction of the dominant 
field distribution. When compared to the conventional linear tetrahedral elements, 
the efficiency of this expansion is apparent (i.e.. many more tetrahedrals are needed 
to model the same cable region). 

6.3.4 Results and Conclusion 


To validate our proposed feed simulation, two circular patch antenna configura- 
tions were used for calculation. One patch antenna was of radius 1.3 cm printed 
atop of a circular cavity (radius=2.1 cm) filled with a dielectric (tr=2.9) material 
0.41 cm deep. For this patch, the feed was placed 0.8 cm from the center and the 
input impedance was measured over the band 2-5 GHz. In fig. 6.6 we compare the 
measured input impedance with data computed on the basis of the proposed equi- 
potential feed model. Clearly, the results from measurements and the equi-potential 
model are in excellent agreement whereas the probe model yields substantially inac- 
curate results near resonance. 

Figure 6.7 shows the comparison between measurements and calculations for an- 
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otlifr [)atrli anleiina wliosr input impcclaiK<‘ was iiicasiircd 1>\ Alxulr am! I’o/ai li > . 

I his patrli liad a radius of 2.0 nil and the O.dlJ' cm thick siilisl latc ha<i ( . = 2.-m and 
a loss tanj^cnt tand = 0.0012. The toed was locati'd O.i cm from tlic (cntci. and U)i 
our FE-Bl calculation the jiatch was placed in a circular cavity of 2.1 1 cm in raiiius. 
As shown in fig. 6.7. the ec|Ui-potential model is again m excellent agia'enu'nt with 
measurements, as opposed to the results by the probe mothd in [6 *]. 

In conclusion, the presented equi-potential feed model has been sliown to lx- 
extremely accurate in modeling coax feed structures. Moreover, its im])lem<'nt at ion 
in the context of a finite element formulation is very simple and as easy to impUnnent 
as the probe feed. It was also demonstrated how the proposed feed model can be 
generalized to the case of asymmetric feed structures where evanescent modes ma\- 
be present. 

6.4 Conclusion 

In developing numerical techniques for antennas, the feed network is one of chal- 
lenging problems to solve in consideration of accuracy, efficiency and simplicity. This 
is primarily because the antenna feed in fabrication has certain instrumental uncer- 
tainties on one side. On the other, as aforementioned, the accuracy of numerical 
results is usually extremely sensitive to the feed model, the feed location, sampling 
rate around the feed point, and so on. 

The proposed numerical feeds in this chapter resemble the practical systems as 
closely as possible, and wdth a thorough consideration of their numerical implemen- 
tations, we realize that they can be used for mostly encountered antenna problems. 
As an addition to the group of feed models, we also developed a circuit modal feed 
wffiich coincides with domain truncations. Since this model has to do with microwav(' 
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Figure 6.3: Illustration of a cavity-backed patch antenna with a coax cable feed. 
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Figure 6.4: (a) Side view of a cavity-backed antenna with a coax cable feed: (b) 
Illustration of the FEM mesh at the cavity-cable junction (the field is set 
to zero at the center conductor surface). 
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Figure 6.5: Field distribution in a shorted coax cable as computed by the finite ele- 
ment method using the expansion (6.18). — : analytical; xxx: numerical, 
(a) Field coefficient cq along the length of the cable (leftmost point is the 
location of the short); (b) Field along the radial coordinate calculated at 
a distance A/4 from the shorted termination. 
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Figure 6.6: Measured and calculated input impedance for a cavity-backed circular 
patch antenna having the following specifications: patch radius r=13mm: 
cavity radius R=21.1mm; substrate thickness t=4.1mm; £^=2.4; and feed 
location x/=0.8 cm distance from center. Results based on the simple 
probe model are also shown for comparison. Our modeling retains the 
vertical wire connection to the patch and uses the incoming coa.xial mode 
field for excitation, (a) Real part; (b) Imaginary part. 
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Figure 6.7: Measured and calculated input impedance for a circular patch an- 
tenna having the following specifications: patch radius r=2cm: substrate 
thickness d=0. 21844cm; feed location from center a/=0.7cm; Cr=2.33; 
fan6=0.0012. [65]. — : measurement; xxx: this method; o o o: probe 
model [65] (a) Real part; (b) Imaginary part. 
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CHAPTER VII 


Circuit Modeling 


Many designs of microstrip antennas require certain feedlines to carry electromag- 
netic signals from the source. Finite element methods are also suited and applicable 
to these wave propagation problems. This chapter is devoted to circuit modeling. 
After an overview of recent mesh termination techniques, we discuss a numerical 
de-embedding process appropriate for the finite element analysis, and then turn to 
the topic of mesh truncations for circuit simulation. 

7.1 Introduction 

One of the most important aspects of finite element implementations is the trun- 
cation of the computational volume. An ideal truncation scheme must ensure that 
outgoing waves are not reflected backwards at the mesh termination surface, i.e. the 
mesh truncation scheme must simulate a surface which actually does not exist. To 
date, a variety of non-reflecting or absorbing boundary conditions (ABCs) have been 
employed for truncating the computational volume at some distance from the radi- 
ating or scattering surface, and applications to microwave circuits and devices have 
also been reported. The ABCs are typically second or higher order boundary condi- 
tions and are applied at the mesh termination surface to truncate the computational 
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voliiuH- as required by an> PDK solution, .\moiiu tlicm. a das> of ,\H( - i- ba-cd un 
tbconr-way \va\c eciuat ion nu'tbod [(•(>, (lij ami am>ibci is dcri\(‘d siaiiins: with tlu' 
Wil ro.\ Lxj)ansiori [GS.dl]. .Also. Inglicr order .\H( ’s using Higdon s Ifi'f. Ttllioi inula- 
tion and problem specific numerical .ABCs have been successfully use'd. jiarl i< iilai h 
for truncating meshes in guided structures [71], There are several difficult i('s with 
traditional .ABCs. .Among them is accuracy control, conformality, ea.se of paralleliza- 
tion and implementation difficulties when dealing with higher order .ABC's. .Also, tlu' 
applications of .ABC's in microwave circuit modeling requires o priori knowh-dge of 
the propagation constants which are typically not available for liigh density packages. 

An alternative to traditional ABCs is to employ an artificial absorber for mesh 
truncation. Basically, instead of an ABC, a thin layer of absorbing material is used 
to truncate the mesh, and the performance for a variety of such absorbers have been 
considered [72, 73]. Nevertheless, these lossy artificial absorbers (homogeneous or 
not) still exhibit a non-zero reflection at incidence angles away from normal. Re- 
cently, though, Berenger [74] introduced a new approach for modeling an artificial 
absorber that is reflectionless at its interface for all incidence angles. In two dimen- 
sions, his approach requires the splitting of the field components involving normal 
(to the boundary) derivatives and assigning to each component a different conduc- 
tivity. In this manner an additional degree of freedom is introduced that is chosen 
to simulate a reflectionless medium at all incidence angles. Provided the medium 
is lossy, this property is maintained for a finite thickness layer. Berenger refers to 
the latter as a perfectly matched layer(PML) and generalization of his idea to three 
dimensions have already been considered [75,76]. Also, implementations of the ab- 
sorber for truncating finite difference-time domain(FDTD) solutions has so far been 
found highly successful. .Nevertheless, it should be noted that Berenger's PML does 
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!iol satisfy Maxwell equations and cannot he easily implemented in linile element 
( FF,M ) solution. 

.-\ new anisotropic(uniaxial ) artificial ahsort)er i] was introduced leu'ntlv foi 
truncating FEM meshes. This artificial absorber is also reflectionless at all incidemc 
angles. Basically, by making appropriate choices for the constitutive ))arameter t<'n- 
sors. the medium impedance can be made independent of frequency, polarization, 
and wave incidence angle. PML layer can then be constructed b\ intioducing suf- 
ficient loss in the material properties. The implementation of this artificial absorber 
for truncating finite element meshes is straightforward and. moreover, the absorber 
is Maxwellian. 


7.2 Numerical De-embedding 


De-embedding presented here is a numerical process used to extract certain circuit 
quantities. Specifically, we are interested in S-parameters for a uniform transmission 
line terminated with any loads denoting the possible discontinuities, which may arise 
from line-to-line or line-to-antenna couplings. The dominant transmission line mode 
is assumed at and near a reference plane Srej in this discussion. 

Consider a transmission line of certain length as shown in fig. 7.1. W ith an 
appropriate shielding scheme, the line is included in the computational domain. The 
full wave analysis provides the E field distribution anyw^here including the region 
along the line. One is therefore able to represent E field along the transmission line 
with respect to the locations to get 

V{z) = Vie~^^ + Vre'^^ (<-l) 

where V' is proportional to the magnitude of E with V, being the incoming and 
the reflected wave amplitude, r is measured from the reference plane Sr,j. y is the 
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Figure 7.1: Illustration of a shielded microstrip line. 


propagation constant to be determined and R = U/\l is the reflection coefficient, 
or 5’n. 

Since in (7.1) 7 . Vj and V'r are three independent quantities that characterize the 
wave propagation and reflection, to determine them one needs to specify the field 
values for V'(c) at three points r_, cq and z+ of equal inter-distances 


V{z.) V{zo) V(c+) 

z+ — zq = zq — z- ( 7 . 2 ) 

along the line. To simplify the problem, we choose the reference plane right at the 
center point such that 20 = 0 and 2 + = — 2 _ = d. Given the three field values from 
FEM computations, it follows 


V{d) = + v;e'"' 

(7.3) 

V''(o) = v: + v; 

(7.4) 

V{-d) = 

(7.5) 


To solve for 7 , we first add (7.3) to (7.5) to get 


(V; + V;)(e^‘' + e-^'') = V-(r/) + V-(-^) 


( 7 . 6 ) 
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from which -• can be determined. The effective guided wavelength A,, and ('fh'ctive 
dielectric constant c^jj may then be calculated by 
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with 3 = and Ao being the free space wavelength. From(7.3) and (7.4). \, 

and V'r are expressed as 

^ V(0)e-^ - V(of) (- 51 

^ 2 sinh( 7 (/) 

v; = V''(o)-v; (f-9) 


Therefore, the reflection coefficient becomes 
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v; 


(7.10) 


This de-embedding process is suited for one port network analysis, ffowever. 
the technique may be readily extended to two-port networks. For instance, on the 
assumption of a perfect termination or match at port 2, once is determined. 52i 
can be obtained by Vo/K , where is the outgoing wave at port 2 predicted by FEM. 

As mentioned before, S-parameter evaluations depend on termination methods. 
Low quality terminations result in prediction errors and make the analysis less reli- 
able. Therefore, high performance termination methods are always desirable and we 


next discuss this issue. 



7.3 Truncation Using DMT 


As alr<‘Hcl\ indicated. S-paianieters from a possiltlc tliscont muit if 2 ,ioii alone, a 
transmission lim* may be e.xtracted at a distant reference' ])lane. where their exists 
only a dominant mode. For shielded microstrip lines at the input i)ort |#1) and 
output port (#2). (similar to that in fig. 7.1). the modes underneath the lines are 


given by 


£'y(-r) = 


Eo{.r)(e + Re' 
r£o(.T)e"^-' 


- e S.„ 


= e Sout 


where Eo{x) denotes the field distribution of incident wave at the incident plane 
(port 1) and R is the reflection coefficient at the same plane. T represents the trans- 
mission coefficient measured at the plane Sout (port 2). and y,. -)2 are the effective 
propagation constants at port 1 and port 2, respectively. 

For truncating the FEM mesh at a specific port, it is necessary to first determine 
the E field pattern across the shielded structure. This can be accomplished by 
assuming a static model shown in fig. 7.2, where the static potential satisfies 

W = 0 


0 = Vo on metallic lii 


(7.12) 


where E = -V0. Sove this standard PDE model, and with a tedious mathematical 


derivation, it is finallv found that 
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A complete FEM system may now be constructed by introducing the ENI fields 
at 5m to truncate the computational domain. This truncation simultaneously in- 
troduces an excitation to the numerical system and the S-parameters may then be 
extracted bv measuring the field distributions at the input and output ports as men- 
tioned before. 

7 A Truncation Using PML 

Below, we begin with a brief presentation of the artificial absorber, and this is 
followed by an examination of the absorber’s performance in terminating guided 
structures and volume meshes in scattering problems. Results are presented which 
show the absorber's performance as a function of thickness/frequency and for differ- 


ent loss factors. 


7.4.1 Theory 
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Figure 7.3: A rectangular waveguide (a) and a microstrip line (b) truncated using 
the perfectly matched uniaxial absorbing layer. 


element method. For a general uniaxial medium, the functional to be minimized is 
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in which /r,. and tr denote the permeability and permittivity tensors whereas E is the 
total electric field in the medium. The surface integrals over 5,n and S'out must be 
evaluated by introducing an independent boundary condition and the .A.BC serves 
for this purpose but alternatively an absorbing layer may be used. An approach to 
evaluate the performance of an absorbing layer for terminating the FE mesh is to 




extract the refiertioii coeffici(‘nt coinj>iitecl in the i)i('>eme ut th<- al>M>rliin;^ l.i\ei U'<'<i 
to terminate the computational domain. Iti thi> stnd\ we (on^id<■l the pet lot nniiu i 
of a thin uniaxial layer for terminating the }■ L mesh in a rectangular uaveguide 
and a microstrijr line. Such a unia.xial layer was projrosed In Sacks il.al. [i ij who 
considered the plane wave reflection from an anisotropic interface (see hg. i . l ). It 



Figure 7.4; Plane wave incidence on an interface between two diagonally anisotropic 
half-spaces. 


J[^ and Cr are the relative constitutive parameter tensors of the form 
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the TE and TM reflection coefficients at the interface (assuming free space as the 
background material) become 


TE 




R™ = 


cosdi — yj ^COsOi 
cos6, + yf^cosOt 


^COsOt — cos6, 
<12 

cosOi + yf^cosOt 
V 02 


7.16) 



and 1)\ clicjosino — U, and r.» = — it follows that ^ = (I foi all iiu td('iMc 
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aiigif's. implying a i)orfertly malchcd niat<Tial in!crfar('. H' u<- mM a, = o - / V. tin- 
rrfiertpd field for a metal-hacked uniaxial layer is 

( 7 . 17 , 

where t is the thickness of the layer and 6, is the plane wave incidence angh'. I lie 
parameter ak is simply the wavenumber in the absorber. Basicallw the pro))os('d 
metal backed uniaxial layer has a reflectivity of -30 dB if .i1ro,-<0, = O.'iT'iA or - 
55dB if 3tcos0, = 0.5A. where A is the wavelength of the background material. 
The reflection coefficient (7.17) can be reduced further by backing the layer with an 
ABC rather than a PEC. However, the PEC backing is more attractive because it 
eliminates altogether the integrals over the surfaces. Clearly, although the interface 
is reflectionless, the finite thickness layer is not and this is also true for Barenger's 
PML absorber. 

Below we present a number of results w'hich show' the performance of the proposed 
uniaxial absorbing layer as a function of the parameter d, the layer thickness i and 
frequency for the guided structures shown in Figure 7.3. \\*e remark that for the 
microstrip line it is necessary to let Ua = fr(,(Q - jp) for the permittivity tensor 
and 02 = - jp) for the permeability tensor, where £^6 and firb are the relative 

constitutive parameters of the background material (i.e. the substrate). 

7.4.2 Results 
Rectangular Waveguide 

Let us first consider the rectangular waveguide showm in fig. 7.3. The guide's cross- 
section has dimensions 4.755 cm x 2.215cm and is chosen to propagate only the T Eu> 
mode. It is excited by an electric probe at the left, and fig. 7.6 shows the mode field 



strciigtli inside the u aveguidc which has Ikmmi terminated 1>> a [uMieetlc matihed 
uniaxial layer. .\s expected, the field ch'cay inside' the ahsurher is exieoiiential ami 
for .i values less than unity the wave does not have suflicic'iit deca\ to sn|)picsv 
reflections from the metal backing of this -bem layer. C’onsec|uent ly. a \ S\\ K of 
about 1.1 is observed for 3 - 0.5. However, as 3 is increased to unity, the \ S\\ H 
is nearly 1.0 and the wave decay is precisely given by ^ 

t is the wave travel distance measured from the absorber interface. F - 2.i//A, and 
here 9, = 44. 5^'. It is noted that when 3 is increased to larger values, the ra]>id 
decay is seen to cause unacceptable VSWR's. One is therefore prompted to look 
for an optimum decay factor for a given absorber thickness and fig. 7.7 provides a 
plot of the TEio mode reflection coefficient as a function of 23tj\g. where we chose 
to normalize with respect to the guided wavelength A^. Figure (.7 is typical of the 
absorber performance and demonstrates its broadband nature and the existence of 
an optimum value of for minimizing the reflection coefficient. Basically, the results 
suggest that 3 must be chosen for a given absorber thickness to provide the slowest 
decay wdthout causing reflections from the absorber backing. That is, the lowest 
reflection may be achieved when the entire absorber width is used to reduce the 
wave amplitude before it reaches the absorber’s backing. As expected, this optimum 
value of 3 changes with frequency but the broadband properties of the absorber are 
still maintained since acceptable low reflections can still be achieved for unoptimized 
3 values. For example, in the case of / = i.5GHz (dashed line) the optimum value of 
3=1 gives a reflection coefficient of -45dB whereas the value of = 3 (corresponding 
to 2/?f/Ag = 2.3) gives a reflection coefficient of -37 dB which is still acceptable for 
many applications. It should be noted though that setting /? = 3 allows use of an 
absorber which is about 2cm or 1/3 free space wavelengths. .Also, as can be realized 
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tlif clisrr('t izfit i(Mi rate plays a roir in fimlins; tlir u|>iinnmi \alu<- uf it A, and iliii' 
the j)resent«‘cl riir\cs rrfer to a sampling rate' of arotiml A. tor thr u a\ ciiti idc 
example. 

Not surprisingl\ (see (G.G)). for this example, the value of o does not play an 
important role in the performance of the absorber and this is demonstrated in fig. 7.>. 
-As seen, setting always a = 3 gives the same performance as the case of n = 1 shown 
in fig. 7.7. Our tests also show that other choices of o give the similar absorlxM 
performance. However, it is expected that a will play a role in the j)res«'me of 
attenuating modes and it is therefore recommended to choose o = J to ensure that 
all modes are absorbed. 

Microstrip line 

The performance of the perfectly matched uniaxial layer in absorbing the shielded 
microstrip line mode is illustrated in fig. 7.9 where the reflection coefficient is plotted 
as a function of 23t/Xg, where \g = Xof^/t^ and Ce// is the effective dielectric 
constant. In this case, the microstrip line is terminated with a 1.87 cm thick. 5 
layered absorber and the line is extended up to 4 layers inside the absorber to avoid 
an electric contact with the metallic wall. Similarly to the waveguide, we again 
observe that an optimum /3 value exists and it was verified that in the absorber 
the wave exhibits the same attenuation behavior as shown in fig. 7.6. The reflection 
coefficient at the optimum /? = 1 is now —4'2dB and if better performance is required, 
a thicker absorbing layer may be required. Again as in the case of the waveguide 
example, the value of q plays little role in the performance of the absorber and this 
is illustrated in fig. 7.10. However, of importance is the behavior of the reflection 
coefficient as a function of 20tlXg. For the waveguide and microstrip examples, we 
observe that the absorption is maximized for approximately the same value of 23t / X,^ 




(al)out O-M. 1 lins ihoc cur\(*s can t»c ux'd fi>r otlici application' a' well, allhouiili 
it should be noted that the discretization rat<' plays an e(piall\ iinpoM.int role and 


tills needs further imestigation. 

The accuracy and validity of the FML applications for circuit iiaranieti'i compu 
tations can also be seen from the result illustrated in fig. 7.11. It is seen that use 
of the optimized 4. -5 cm PML layer, with o = 1 and .f = 1. yields ver\ accuiate 
input impedance values. The shown microstrip line impedances were computed by 
measuring the vertical field at the probes location without a need to extiact tlu 
VSWR which is often difficult with unstructured finite clement meshes. Note that 
the shielded microstrip line dimensions for the data are given in fig. i.ll. 
Meanderline 

Another example is the meander line shown in fig. 7.12. For the FEM simulation, 
the structure was placed in a rectangular cavity of size b.Srnm x 18.0mm x 3.175mm. 
The cavity was tessellated using 29 x 150 x 5 edges and only 150 edges were used 
along the y-axis. The domain was terminated with a 10 layer PML. each layer being 
of thickness i = 0.12mm. The Sn results are shown in fig 7.13 and are in good 
agreement with the measured data [78]. 
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Figure 7.5; Typical field values of TEio mode inside a rectangular waveguide termi- 
nated by a perfectly matched uniaxial layer. 



Figure 7.6: Field values of the TEio mode inside a waveguide terminated by a per- 
fectly matched uniaxial layer. The absorber is 10 elements thick and each 
element was 0.5 cm which translates to about 1.3 samples per wavelength 
at 4.5 GHz. 












Figure 7.9: Reflection coefficient vs 2(31 jXg with q= 1, for the shielded microstrip line 
terminated by the perfectly matched uniaxial layer. 



Figure 7.10: Reflection coefficient vs 2j3tfXg with q = for the shielded microstrip 
line terminated by the perfectly matched uniaxial layer. 
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Figure 7.11: Input impedance calculations for the PML terminated microstrip as 
compared to the theoretical reference data. 











CHAPTER VIII 


AWE: Asymptotic Waveform Evaluation 
8.1 Brief Overview of AWE 

Although full wave electromagnetic systems are large and cumbersome to solve, 
typically only a few parameters are needed by the designer or analyst. A reduced or- 
der modeling of these parameters (input impedance, S parameters, far field pattern, 
etc) is therefore an important consideration in minimizing the CPU requirements 
needed for generating the frequency response of the parameter. The Asymptotic 
Waveform Evaluation( AWE) method is one approach to construct a reduced order- 
ing model of the input impedance or other useful electromagnetic parameters. .'\WE 
relies on a Fade approximation of the given parameters to avoid the repeated solution 
of the system at each frequency value. It has already been applied to problems in cir- 
cuit analysis and in this paper we demonstrate its application and validity when used 
in conjunction with the finite element method to simulate full wave electromagnetic 
problems. 

The method of Asymptotic Waveform Evaluation (AWE) is a reduced-order mod- 
eling of a linear system and has already been successfully used in VLSI and cir- 
cuit analysis to approximate the transfer function associated with a given set of 
ports/ variables in circuit networks [79-82]. The basic idea of the method is to de- 
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ail aj)j)roxiniat e iraiisfer function of a cix'cn linear >y>iem fixini a liiniteil >ei i>f 
spectral solutions. Topically, a Fade expansion of the trati^fei function i> posiiilaled 
whose coefficients are then determined i»\ mat chine, the Fade rejire-'ent at ion to the 
available spectral solutions of the complete system. 

In this chapter we investigate the applicability of the .-\\\'E method for approx- 
imating the response of a given parameter in full wave simulation of radiation or 
scattering problems in electromagnetics. Of particular interest is to use .\\\'E for 
evaluating the input impedance of the antenna over an entire bandwidth from a 
knowledge of the full wave solution at a few (even a single) frequency points. .Also, 
the method can be used to fill-in a backscattering pattern with respect to frequency 
from a priori knowledge of the simulation system with a few data samples of that pat- 
tern. Given that practical partial differential equation(PDE) systems involve several 
thousand unknowns, .AWE can indeed have a dramatic reduction of CFU require- 
ments in generating a response for a given system parameter (state variable) without 
a need to resolve the system for the fields in the entire computational grid. Below we 
first describe the recasting of the FEM system for application of the AWE. We then 
proceed to describe the AWE method and demonstrate its application, accuracy and 
efficiency for computing the input impedance of a shielded microstrip stub. 

8.2 Theory 

8.2.1 FEM System Recast 

The application of the finite element method to full wave electromagnetic solu- 
tions amounts to generating a linear system of equations by extremizing the func- 
tional [83] 
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\vlicr<’ <•> denote^ an iiiiK'i-product and b.t. i> a ])<.)'sit>lc l)i>iuidar\ t<'iin v.lm'c 
specific form is not recjuircd for tins discussion. .-\lso. llic dyadics o and ■>' arc mate- 
rial related coefficients. A- = — = — is tlie waveniuiiher and is llu- ( i)rresi)oiuiin!i 

A c 

ojrerating frequency with c being the speed of light. .-\ disen-tized form of (S.li 
incorporating the appropriate boundary conditions is [29] 

(Ao-fA-A, +eA,]){.V) = {/) (sa 

wliere A, denote the usual square (sparse) matrices and {/} is a column matrix 
describing the specific e.xcitation. 

Clearlv (8.2) can be solved using direct or iterative methods for a given value of 
the wavenumber. Even though A, is sparse, the solution of the system (8.2) is com- 
putationally intensive and must be further repeated for each k to obtain a frequency 
response. Also, certain analyses and designs may require both temporal and fre- 
quency responses placing additional computational burdens and a repeated solution 
of (8.2) is not an efficient approach in generating these responses. .An application of 
AWE to achieve an approximation to these responses is an attractive alternative and 
below we formulate AWE in connection with the FEM system (8.2), whose imple- 
mentation is considered in connection with antenna and microwave circuit problem. 
For these problems it turns out that the excitation column {/} is a linear function 
of the wavenumber and can therefore be stated as 

{/}=*{/,) (8.3) 

with {/i) being independent of frequency. This observation will be specifically used 
in our subsequent presentation. 
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8.2.2 Asymptotic Waveform Evaluation 

To dcscrilu- tlic basic idea of AWE in conjunction willi the 1 I'.M. we Ix'iiin b\ 
first expanding tlie solution {A'} in a laylor series alunit k\, a^ 

{ A } = {Ao} + (A- — A'o) {Ai } + (A- — Ao)‘ {A j} + . . . 

+ (^‘ ~ A'o)^ ~ Ao)^’*"' } (S. 1 ) 

where {A'o} is the solution of (S.2) corresponding to the wavenutnlx'r Ay. By intro- 
ducing this expansion into (8.2) and equating equal powers of A- in conjunction with 
(8.3). after some manipulations we find that 

{A-„) = <.-oA„-'{/,} 

{A',} = A„-‘l{/,)-A,{.V„}-2fc„AH.Vo}] 

{-^ 2 } = ~Ao {-^ 1 } + -A.2({Ao} + 2A:o {Ai})] (8.5) 


{A/} — — Aq ^ [Ai {A/_i } + A2({A;_2} + 2A:o {A/_i ) )] 

with 

Ao = Ao + AtoAi + A;qA2 (S-6) 

Expressions (8.5) are referred to as the system moments whereas (8.6) is the 
system at the prescribed wavenumber (A:o). Although an explicit inversion of Aq ' 
may be needed as indicated in (8.5), this inv^ersion is used repeatedly and can thus 
be stored out-of-core for the implementation of .AWE. Also, given that for input 
impedance computations we are typically interested in the field value at one location 
of the computational domain, only a single entry of {A'/(A’)} needs be considered, say 
(the pth entry) Xf(k). The above moments can then be reduced to scalar form and 



t 1 j(‘ rx])cinsions ( ) lK*roino a scalar r('prcstMil at ion of A < ) al>out t hr < on hum iiu 

solulioii at Aq. Io yield a more cDiivergent exj)r(*ssioii. \vc‘ can in>trad vc\cv\ ihr 
moments to a Fade expansion, whicli is a coinenlioiial rational iunction in form. A 
special case of the q\U order of such an expansion is given 1)\ 


rp I ^^0 + “ ^'o) + ^^2(A' — Aq)^ -h . . . + A — AoF 

^ ^ 1 + 6]( A Ao) + i> 2 (A — Ay)*-’ + ... + 6.^(A — A^)^ 


(S.7) 


where czj and 6^ {? — 0. 1 q) are referred to as the Fade coefficients. 

For transient analysis, it is observed that the Fade expansion can be rc'formulat (h 1 
by partial fraction decomposition [82,84] as 


= A-j, + 


1 = 1 


k — ko — k, 


(S.8) 


where A',o >s the limiting value when k tends to infinity. Clearly, this is the represen- 
tation suitable for time/frequency domain transformation. The residues and poles 
(r, and ko + C) in (8.7) or (8.8) correspond to those of the original physical system 
and play important roles in determining the accuracy of the approximation. In gen- 
eral, higher order expansion contains more system residues and poles and usually 
provides a better approximation. Since the accuracy of AWE relies on the dominant 
residues and poles located in a complex plane closest to the point on the real axis 


ko from the origin, in practice the number of poles (and residues) needed to obtain 
a sufficiently accurate expansion can be much smaller than that of the original nu- 
merical system, which is the beauty of AWE method. (Detailed analysis and theory 
of Fade expansion can be found for instance in [85].) 

For hybrid finite element - boundary integral system, the implementation of 
AWE is more involved because the fully populated submatrix of the overall system 
may be associated with a more complex dependence on frequency. In this case it is 
attractive to instead generate the full submatrix by introducing a spectral expansion 



of ihe ox[>oiKMit ia! Wouiidary intenral kt'riu‘l lu faciliiair t hr r\i rad ion <>1 t 1m* 
niomenls. J hi> a[)proach does increase ilu' ('omplic ai ion> for impiriiK'in iim AWl ., 
It however remains far more efficient in term> of ( IH laxpiiiannent^ wIk'h ('ompan'd 
to the conventional approach to continuously rep(xuing tin* solution of tlu' eniii(‘ 
system. 

8.3 Numerical Implementation 

As an application of AWE to a full wave electromagnetic simulation, we consider 
the evaluation of the input impedance for microstrip stub shielded in a metallic 
rectangular cavity as shown in fig. 8.1, As expected, the stub's input impedance is 
a strong function of frequency from 1-3 GHz and this example is therefore a good 
demonstration of AWE’s capability. 

The shielded cavity is 2.38cm x 6,00cm x 1.06cm in size and the microstrip stub 
resides on a 0.35cm thick substrate having a dielectric constant of 3.2. The stub is 
0.79cm wide and A/2 long at about 1.8 GHz. We note that the cavity is terminated 
at the perfectly electric conductor (PEC) back wall by an artificial absorber having 
relative constants of Cr = (3,2, —3.2) and fir = (1.0.— 1.0). In this study the artificial 
absorber was used for setting up an appropriate forced problem rather than to es* 
tablish a perfectly matched interface. Nevertheless the numerical FEM system was 
already demonstrated valid and accurate for microwave circuit analysis [86]. 

The frequency response of the shielded stub was first computed using a full wave 
finite element code from 1 to 3 GHz at 40MHz intervals (50 points) to serve as the 
reference solution. We then chose the single input impedance solution at 1.78GHz 
in conjunction with the 4th order and 8th order AWE representation given in (8. 8) 
to approximate the reference response. As seen in fig. 8.3. the 4tli order AWE 
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roprcx-iitat ion in aurccrnrnl witli tlir r(‘al and reactive part- i>l lln' ifteiciK e inpni 
iiiip(‘daii(‘(‘ solnl ion o\ei about )(>' ( and ddV< baiuKvult !i. ropect i\ cl\ . 1 hi> ileail\ 

slious that the contributions of tlie systeni poles in the complex /r |>lane lead lu an 
accuracy difference to the real and reactive coni])onents. Surprisinp.l\ . the mIi oidt'i 
AW E representation recovers the reference solution over the entire 1 - KiH/ band for 
both impedance components. We also observed that the ( PI recpnreim'iil s foi 1th 
and 8th order computations are nearly the same except for a few more tinu's of 
matrix-vector products. The number of these product ojierations is in the order of 
the AWE approximation order q and therefore much smaller than the size of the 
original numerical system. 

It is also apparent that to demonstrate the AW'E efficiency we only solved the 
system once at one frequency point. The save of CPU time can be easily estimated 
when compared to solve the system conventionally for each frequency over the entire 
band. Thus, the AWE representation is an extremely useful addition to electro- 
magnetic simulation codes and packages when a wideband frequency response of the 
system is required. The development and utility of the method for more complex 
numerical systems and multiple parameter simulation can be readily extended and 
will be considered in the future. 
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Figure 8.1: Illustration of the shielded microstrip stub excited with a cunent probe. 



Figure 8.2: Impedance calculations using traditional FEM frequency analysis for a 
shielded microstrip stub shown in figure 8.1. Solid line is the real part 
and the dashed line denotes the imaginary part of the solutions. These 
computations are used as reference for comparisons. 
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Figure 8.3: 4th order and 8th order AWE implementations using one point expan- 
sion at 1.78 GHz are shown to compare with the reference data. With 
the 4th order AWE solutions, 56% and 33% bandwidth agreement can be 
achieved for the real (a) and imaginary (b) parts of impedance compu- 
tations, respectively. It is also shown that the 8th order solutions agree 
excellently with the reference data over the entire band, (a) Real Part 
(b) Imaginary Part computations 








CHAPTER IX 


Conclusions 

9.1 Discussion on the Research Work 

During the period of developing the hybrid finite element methods, many ex- 
pected and unexpected issues were frequently encountered. .Among them are the 
understanding of physical systems, development of mathematical models, interpre- 
tation of results, lack of measurement data for comparison, and increased compu- 
tational demands, etc. VVe can comfortably state that significant progress has been 
made during the course of this work. Some of our accomplishments are summarized 
below. 

• General purpose hybrid FE-BI method development 
Once the FE and BI subsystems and the hybrid method were mathemati- 
cally formulated, a major effort was then devoted to the integration of the 
two subsystems. The interface between the FE-BI program and a commercial 
(SDRC-IDEAS) mesh generator was developed with minimum but sufficient 
geometry and meshing data. The latter task was important in permitting the 
geometrical modeling and meshing of printed antenna configurations of arbi- 
trary shape. It is this general version of the FE-BI code that can (in theory) 
be used to simulate any planar conformal antenna. 
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• ITI;KAT1\ t SYSTEM SOIAT.H 

A mriiiorN’ saving algorithm II PACK was intrrt wiiiccl witli tlic li\lni»l 1 IM 
sul)syslem to registrr only the non zero FKM entries. 1 he Hi( (> iteiatiw 
solver was independent Iv de\eloped for partially s))arse and paitiall\ full ma- 
trices in conjunction with the ITPAC K algorithm. 

• Uniform grid BI subsystem — BiCG-FFT 

To facilitate the efficient storage and evaluation of the BI sub-system, a uni- 
form right triangular zoning scheme for discretization of the boundary integral 
equation was introduced by re-numbering the triangle edges as dictated by their 
geometrical locations. This approach leads to a BI sub-system which could be 
cast as a 2-D discrete convolution, thus allowing use of FFT for fast execution 
of the iterative solver. This truncation/termination the "exact" evaluation of 
rectangular and right-triangular patches. 

• Non-uniform BI subsystem — Overlay-BiCG-FFT 

For non-rectangular patches, an interpolation scheme was proposed to make use 
of the efficient BiCG-FFT technique by overlaying a fictitious uniform grid with 
the original arbitrary mesh. The forward/backward transformation matrices 
to account for field interpolations using localized basis functions were derived 
and they were indeed highly sparse. 

• Feed modeling 

Feed modeling is one of the most important and challenging tasks in the context 
of the general purpose FEM. To this end, a series of commonly used feed 
structures were modeled using the hybrid technique, especially in consideration 
of efficiency and accuracy. These include probes/generators, aperture coupled 
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slotliiic. inicroslrip Wuv. (oax cable, etc. 

• Prismatic FEM Ellmcnts Incorporation 

A major problem in an\ hybrid FEM analysis is \]\v t(‘diou> prr-piucosiim 
for mesh generation. Thin layer substrates in the j)resenc(‘ of thick s[)a((M(sl 
are often found in practical conformal antenna designs. Howev^u. this typical 
configuration leads to large numerical systems when tetrahedral elements are 
used. To alleviate these difficulties, the prismatic edge-based elements were* 
developed and incorporated in the hybrid system. This formulation exhibits 
certain features/advantages that tetrahedral FEM does not. It can therefore 
be used to compensate the tetrahedral FEM as a subsystem module. 

• Mesh Truncations With DMT and PML 

The uniaxial or other anisotropic medium simulation may be readily accom- 
plished using the proposed hybrid FEM technique due to the geometrical adapt- 
ability of the tetrahedral elements. Hence the PML was first introduced into 
the 3‘D FEM. Various performance studies were carried out to optimize the 
application of the PML to microwave circuit simulations. In the meantime, an 
analytical approach, dominant mode truncation (DMT), was proposed and im- 
plemented as an alternative mesh truncation of the FEM domain for microstrip 
lines and similar structures. 

• Reduced order approximation — AWE 

AWE has been reported useful in RLC and VLSI applications. For wideband 
and highly varying frequency responses, this technique is particularly efficient. 
Given the promise of the method for broadband simulations of VLSI circuits, 
we consider its application to electromagnetic system. In particular. AWE 



was iiirorporatcd into tlic finite element method. It was indeed oh-eiAcd that 
the attracti\e features of AW E are maintained wIkmi ns<'d in eleii nuiiagnet ic 
problems. 

9.2 Suggestions for Future Tasks 

The following is a list of suggested tasks for further development of the finite 
element methods 

• Higher order edge-b.^sed FEM development 

• Adaptive elements 

• Mixed Elements and interface 

• Anisotropy (with loss) FEM Investigations/applications 

• Incorporation of More Robust truncations 

• Modular development and integration (with user interface) 

9.3 Modular Development 

Hybrid finite element methods for the analysis of various electromagnetic prob- 
lems encountered in practice are still on the way to reach its maturity. As well known, 
any general purpose technique (such as the commercially available software in electro- 
magnetics) either loses its efficiency or becomes incapable when simulating intricate 
problems. It is anticipated that at the current stage of the FEM development with 
the limited capacity of computing resources, more and more specialized techniques 
wdll be desired, particularly when efficiency and speed become a key consideration 
in large scale computations and in engineering design. 



1 1 (. 


With a wiiolc sot of specially developed tecliiii(|ue> ami met liudoluuic'. one '•lioiikl 
then consider to create an integration eiiviroiiim'nt . As shown in lig. 't.l. we pidpo^e 
this FEM modular en\ ironment for future computational elect romagiuM i( a|)pliea- 
lions. A well designed modular finite element methods will he the most (ai)al)l(' and 
robust in the future! 


FEM Multi-Module Environment 



Figure 9.1: Multi-modular FEM environment 
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APPENDIX A 


Evaluation of Matrix Elements for Tetrahedrals 


Referring to Fig. A.l and the associated table, the fields in the (th tetrahedron 
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Table 
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Figure A.l: (a) .A tetrahedron, (b) its local node/edge numbering scheme 
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i=l 

where the basis functions W' are given l>y 

j f:-. + g?-, ref; 

I 0 outside element 

— r, X r,j r,,. r,j : f)osition vectors of vertices /i and i < (sei' lab 

6\^ 

6,fer-. 


(r.2 - r.. ) 

b, 

= |r,j - r,j I = length of the ith edge (see Table) 

= element's volume 

We note that 

V . = 0 y X 'w^ = 2g. 

indicating that W- are divergenceless. Furthermore, 

( 1 1 = 3 
[ 0 f / j 

where r-' has its tip on the jth edge of the tetrahedron. This last property ensures 
that the coefficients = E • e, represent the average field value at the ith edge of 
the tetrahedron. 

Using the above basis functions, we now proceed with the derivation of the matri.x 
elements We have 

[ff l(v X w') ■ (V X w') = -g, ■ gy ; 

J J Jv. 


Wl_.(r) 
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g 7 -, 

e, 

b, 

v; 
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Also. 
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where 

D = (f. X gj) + (fj g,) 

and 


h 

I 2 



f. f.dr 

r ■ D dv 

(g. X r) ■ i&j X r)dv 


Since f is a constant vector. 7i reduces to 


/l = f. ■ f, V e 


To evaluate I 2 we first set 

444 
X = ^ T.x,; y = ~ 

1=1 1=1 t=i 

where x^, y,, [i = 1,... ,4) denote the {x,y.z) coordinates of the tetrahedron s 
vertices and L, are the simplex coordinates or shape functions for the same element. 
That is, L, is the normalized volume of the tetrahedron formed by its three corners 
other than the zth, and the point (x,y,z) located wdthin the tetrahedron. Using the 
standard formula for volume integration within a tetrahedral element and simplifying. 


/2 = 


Vv 


Dx ^ X, + Dy ^ y, 4- D; ^ z, 


1=1 


1=1 


1=) j 


we have 



\\]\orc Dry, is t he nM h roin])OiKMil of D. *1 1h‘ oxaluat iuii ut /. can Ik- >iiu|)lilir(i 1>\ ilu- 
uso of basic vector identities. \\(‘ lia\(' 


h = g, g, ^ (g, r)(g^ -rlr/r 

dv + [Qixdj-r + .9iri/jr) j </ + [fliTfljj “t" 

^ xydv - (9,r</jc + ^ </'■ - ^ 

where represents the 77)th component of the vector g,. Eacli of tlie above integ,i als 
can be easily evaluated analytically and the result can be expressed in the g('uei-,d 
form 


= [g<y9jy ^ 9>--9jz] ^ 

~ d" 9jr9^y) 
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aiUm ^ ~- 


20 


I ^ ^ “b ^ ^ ^ ^ ^t7 

Ltrrl 2 = 1 1 = 1 

for l,m ^ K... .3. The parameters ai or can represent any of the rectilinear 


variables x, y, z. 



APPENDIX B 


Evaluation of the Boundary Integral System 

Matrix 


The explicit representation of the boundary integral subsystem matrix is gi\en 
by (3.12) and can be rewritten as 

Bpq^^ = 2 {-kls, • S, + V X S. ■ V X Sj) Go(r. r') dS' dS (B.l ) 

where Go(r,r') is the free space Green’s function, and Tf is the pth triangle of the 
triangle pair Sp as shown in fig. 3.3. Similar to the finite element assembly procedure, 
it should be recognized that the definition of (B.l) virtually involves an assembling 
over the triangles. 

To proceed, (3.14) is used to discritize the field region and thus its curl is given 

by 


V X Si(r) = t(r)-^£ 


:b.2) 


where e(r) is defined by (3.15). Note that when deriving (B.2), the fact that r is 
located inside the pth triangle in a planar surface is considered and therefore V r = 2. 
Given the Green’s function, it is straightforward to express the matrix entries as 
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dS' dS 
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^JIJL 
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:b.3) 



in which H = jr ~ r^. 1 hcse^ inl(‘"ral> can Ix' r(‘adil\' t^xaluatcd lor non ''clt icll 

terms }>\‘ minierical iiiK'grations. ll is also ol)s('r\(‘d that onct' 7^ coincidi'^ witii 
Tj . the integrands become singular l>ecause of tlu‘ Cireen s fuiution. In this cas(\ 
the singularity should be removed. For the second integral, this is ac(H)mplished b\ 
subtracting and adding an additional term. That is 

ff ff — \ 

— Mr).,(r)<i5V5 = rf>' 

The first integral in (B.4) is evaluated using numerical integrations and the second 
one is carried out analytically [55]. Similarly, the first integral in (B.3) is rearranged 



/ L I Ly 

+ JJ ^(r-r^]-(r-r'^Mr)c,(r)^dS’<IS (B.5) 

in which the first integral on the right hand side is numerically integratable with 
singularity removed and the second one again may be expressed in a simple anal\ tical 


form [55]. 



APPENDIX C 


Formulation for Right Angle Prisms 


For FEM implemenlation. the following quantities are required 


P" = 

mn 

[ V V xV^clV 

(CM) 
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^ v„, • V„ dV 

(C.2) 


where the curls are given by 
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^ ^ + + + Ar-r)] j = 4.5.6 

V X V*; = [(x^2 — 3:;:i )i + (t/fc2 — Vki )y] k = 7,8.9 

To this end, we follow the notation defined in (4.13) and (4.14). where ?. ?'=1.‘2.3 
represent the top triangle edges, j,/=4,5,6 denote the bottom triangle edges and 
k,k'-7,8,9 stand for the vertical three edges. It is found that (C.2) and (C.3) can 
be analytically evaluated and we tabulate the results as follows 
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The remaining quantities in the above list of the expressions are defined as 
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These integrals can be expressed in terms of tlie global coordinates of the three nodes 
(-V,, y,). (A'j. ). (A'n,. Vm )• Specifically, assuming that the three nodes i.j and in of 

a triangle are in counterclockwise rotation, we then have. 

1 y, 

1 Vj 

f it/m 

= J^ xdxdy = Y ) 

S}' = I ydxdy = ^ (V' + V' + V;„) 
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APPENDIX D 


System Derivation From A Functional 


Referring to fig. 2.2. we begin with the functional 

J’(H) = - [ fv X H ■ X H - A-qH • ■ h) dn 

*" J J -\-Ha 


(D.l) 


to derive the system in terms of the scattered field. On inserting the field decompo- 
sition H = the functional becomes 


+ (D.2) 

where the fact that does not exist in has been considered, and (, ) represents 

the integral of the same form as in (D.l). Once a self-adjoint system operator is 
assumed, it then follows that 


( H.nc ^ ) I ) 1 

Also, in free space D/, 

(H-c 

Upon invoking the divergence theory, we have 
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Kvidciit ly. for a s(‘lf*a(ljoinl o[XTatur. onr r<‘a(iily n*( ! hr >\ "Wmii < 2 . I J • rhi <iiiu d 

via (lalerkin s inelhod. It should Ix‘ nolrd that Ix'sidr^^ t lir huundar\ t r<ui>it ion 


conditions, the self-adjoint property of a system operatoi >imply rr(|uir(‘> 



In the case of a noii-self-adjoinl operator, it is generally not |)ossihl(‘ to Kxover tin' 
system given by (2.42) in the same manner. This is l)ecause the functional (D.2) in 
terms of the scattered field is of the form 

JT(H) = i / (v X - 1;' • V X • Ji, ■ dil 

+ \ J (v X • V X H'"" - A-^H""“' • ■ H'"'-) f/Q 

^ \ j (v X H’"'^ • • V X • H*--'"') t/n 


+ / • 

(^hi X 

• V 

X 

1 dS 


Jr, 

(^2 X 

• V 

X 

1 dS 

(D.7 


It is observed that the first integral shows the same form of the FEM system as that in 
(2.42). All other integrals in (D.7) contribute to the system excitation. .Apparently, 
the two integrals over domain are not identical, leading to a different FEM system 
than (2.42). 
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